Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 2007-2009, 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 : /* $Id$ */
16 :
17 : //-------------------------------------------------------------------------
18 : // Implementation of the on-the-fly v0 finder for ITS tracker MI
19 : // Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
20 : // Extraction from AliITStrackerMI: Andrea Dainese
21 : // Current support and development:
22 : //-------------------------------------------------------------------------
23 :
24 : #include <TMatrixD.h>
25 : #include <TTree.h>
26 : #include <TString.h>
27 : #include <TRandom.h>
28 : #include <TTreeStream.h>
29 :
30 :
31 : #include "AliLog.h"
32 : #include "AliESDVertex.h"
33 : #include "AliESDEvent.h"
34 : #include "AliESDtrack.h"
35 : #include "AliESDV0Params.h"
36 : #include "AliV0.h"
37 : #include "AliHelix.h"
38 : #include "AliITSRecPoint.h"
39 : #include "AliITSReconstructor.h"
40 : #include "AliITStrackerMI.h"
41 : #include "AliITSV0Finder.h"
42 : #include "AliKFParticle.h"
43 : #include "AliKFVertex.h"
44 :
45 118 : ClassImp(AliITSV0Finder)
46 :
47 : Bool_t AliITSV0Finder::fgDisabled = kFALSE;
48 :
49 0 : AliITSV0Finder::AliITSV0Finder():
50 0 : fDebugStreamer(0)
51 0 : {
52 : //Default constructor
53 :
54 0 : if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0)
55 0 : fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
56 0 : }
57 :
58 : //------------------------------------------------------------------------
59 : void AliITSV0Finder::UpdateTPCV0(const AliESDEvent *event,
60 : AliITStrackerMI *tracker)
61 : {
62 : //
63 : //try to update, or reject TPC V0s
64 : //
65 0 : if (fgDisabled) {
66 0 : AliWarningClass("On-fly V0 finded is disabled");
67 0 : return;
68 : }
69 0 : const Float_t kMinTgl2= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinTgl2();
70 :
71 0 : TObjArray *trackHypothesys = tracker->GetTrackHypothesys();
72 :
73 0 : Int_t nv0s = event->GetNumberOfV0s();
74 0 : Int_t nitstracks = trackHypothesys->GetEntriesFast();
75 :
76 0 : for (Int_t i=0;i<nv0s;i++){
77 0 : AliESDv0 * vertex = event->GetV0(i);
78 0 : Int_t ip = vertex->GetIndex(0);
79 0 : Int_t im = vertex->GetIndex(1);
80 : //
81 0 : TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)trackHypothesys->At(ip):0;
82 0 : TObjArray * arraym = (im<nitstracks) ? (TObjArray*)trackHypothesys->At(im):0;
83 0 : AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
84 0 : AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
85 : //
86 : //
87 0 : if (trackp){
88 0 : if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
89 0 : if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
90 0 : if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
91 : }
92 : }
93 :
94 0 : if (trackm){
95 0 : if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
96 0 : if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
97 0 : if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
98 : }
99 : }
100 0 : if (vertex->GetStatus()==-100) continue;
101 : //
102 0 : Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
103 0 : Int_t clayer = tracker->GetNearestLayer(xrp); //I.B.
104 0 : vertex->SetNBefore(clayer); //
105 0 : vertex->SetChi2Before(9*clayer); //
106 0 : vertex->SetNAfter(6-clayer); //
107 0 : vertex->SetChi2After(0); //
108 : //
109 0 : if (clayer >1 ){ // calculate chi2 before vertex
110 : Float_t chi2p = 0, chi2m=0;
111 : //
112 0 : if (trackp){
113 0 : for (Int_t ilayer=0;ilayer<clayer;ilayer++){
114 0 : if (trackp->GetClIndex(ilayer)>=0){
115 0 : chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
116 0 : trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
117 0 : }
118 : else{
119 0 : chi2p+=9;
120 : }
121 : }
122 0 : }else{
123 : chi2p = 9*clayer;
124 : }
125 : //
126 0 : if (trackm){
127 0 : for (Int_t ilayer=0;ilayer<clayer;ilayer++){
128 0 : if (trackm->GetClIndex(ilayer)>=0){
129 0 : chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
130 0 : trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
131 0 : }
132 : else{
133 0 : chi2m+=9;
134 : }
135 : }
136 0 : }else{
137 : chi2m = 9*clayer;
138 : }
139 0 : vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
140 0 : if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10); // track exist before vertex
141 0 : }
142 :
143 0 : if (clayer < 5 ){ // calculate chi2 after vertex
144 : Float_t chi2p = 0, chi2m=0;
145 : //
146 0 : if (trackp&&TMath::Abs(trackp->GetTgl())<kMinTgl2){
147 0 : for (Int_t ilayer=clayer;ilayer<6;ilayer++){
148 0 : if (trackp->GetClIndex(ilayer)>=0){
149 0 : chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
150 0 : trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
151 0 : }
152 : else{
153 0 : chi2p+=9;
154 : }
155 : }
156 0 : }else{
157 : chi2p = 0;
158 : }
159 : //
160 0 : if (trackm&&TMath::Abs(trackm->GetTgl())<kMinTgl2){
161 0 : for (Int_t ilayer=clayer;ilayer<6;ilayer++){
162 0 : if (trackm->GetClIndex(ilayer)>=0){
163 0 : chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
164 0 : trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
165 0 : }
166 : else{
167 0 : chi2m+=9;
168 : }
169 : }
170 0 : }else{
171 : chi2m = 0;
172 : }
173 0 : vertex->SetChi2After(TMath::Max(chi2p,chi2m));
174 0 : if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20); // track not found in ITS
175 0 : }
176 0 : }
177 : //
178 0 : }
179 : //------------------------------------------------------------------------
180 : void AliITSV0Finder::FindV02(AliESDEvent *event,
181 : AliITStrackerMI *tracker)
182 : {
183 : //
184 : // V0 finder
185 : //
186 : // Cuts on DCA - R dependent
187 : // max distance DCA between 2 tracks cut
188 : // maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
189 : //
190 : const Bool_t kCheckPropagate = kFALSE;
191 16 : const Float_t kMaxDist0 = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxDist0();
192 8 : const Float_t kMaxDist1 = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxDist1();
193 8 : const Float_t kMaxDist = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxDist();
194 8 : const Float_t kMinPointAngle = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinPointAngle();
195 8 : const Float_t kMinPointAngle2 = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinPointAngle2();
196 8 : const Float_t kMinR = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinR();
197 8 : const Float_t kMaxR = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxR();
198 8 : const Float_t kMinPABestConst = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinPABestConst();
199 8 : const Float_t kMaxRBestConst = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxRBestConst();
200 8 : const Float_t kCausality0Cut = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetCausality0Cut();
201 8 : const Float_t kLikelihood01Cut = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetLikelihood01Cut();
202 8 : const Float_t kLikelihood1Cut = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetLikelihood1Cut();
203 8 : const Float_t kCombinedCut = AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetCombinedCut();
204 8 : const Float_t kMinClFullTrk= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinClFullTrk();
205 8 : const Float_t kMinTgl0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinTgl0();
206 8 : const Float_t kMinRTgl0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinRTgl0();
207 :
208 8 : const Float_t kMinNormDistForbTgl0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinNormDistForbTgl0();
209 8 : const Float_t kMinClForb0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinClForb0();
210 8 : const Float_t kMinNormDistForb1= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinNormDistForb1();
211 8 : const Float_t kMinNormDistForb2= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinNormDistForb2();
212 8 : const Float_t kMinNormDistForb3= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinNormDistForb3();
213 8 : const Float_t kMinNormDistForb4= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinNormDistForb4();
214 8 : const Float_t kMinNormDistForb5= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinNormDistForb5();
215 8 : const Float_t kMinNormDistForbProt= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinNormDistForbProt();
216 8 : const Float_t kMaxPidProbPionForb= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxPidProbPionForb();
217 :
218 8 : const Float_t kMinRTPCdensity= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinRTPCdensity();
219 8 : const Float_t kMaxRTPCdensity0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxRTPCdensity0();
220 8 : const Float_t kMaxRTPCdensity10= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxRTPCdensity10();
221 8 : const Float_t kMaxRTPCdensity20= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxRTPCdensity20();
222 8 : const Float_t kMaxRTPCdensity30= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxRTPCdensity30();
223 :
224 :
225 8 : const Float_t kMinTPCdensity= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinTPCdensity();
226 8 : const Float_t kMinTgl1= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinTgl1();
227 :
228 8 : const Float_t kMinchi2before0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinchi2before0();
229 8 : const Float_t kMinchi2before1= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinchi2before1();
230 8 : const Float_t kMinchi2after0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinchi2after0();
231 8 : const Float_t kMinchi2after1= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMinchi2after1();
232 8 : const Float_t kAddchi2SharedCl= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetAddchi2SharedCl();
233 8 : const Float_t kAddchi2NegCl0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetAddchi2NegCl0();
234 8 : const Float_t kAddchi2NegCl1= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetAddchi2NegCl1();
235 :
236 8 : const Float_t kSigp0Par0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetSigp0Par0();
237 8 : const Float_t kSigp0Par1= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetSigp0Par1();
238 8 : const Float_t kSigp0Par2= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetSigp0Par2();
239 :
240 8 : const Float_t kSigpPar0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetSigpPar0();
241 8 : const Float_t kSigpPar1= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetSigpPar1();
242 8 : const Float_t kSigpPar2= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetSigpPar2();
243 8 : const Float_t kMaxDcaLh0= AliITSReconstructor::GetRecoParam()->GetESDV0Params()->GetMaxDcaLh0();
244 :
245 8 : if (fgDisabled) {
246 0 : AliWarningClass("On-fly V0 finded is disabled");
247 0 : return;
248 : }
249 :
250 :
251 8 : TObjArray *trackHypothesys = tracker->GetTrackHypothesys();
252 8 : TObjArray *bestHypothesys = tracker->GetBestHypothesys();
253 8 : TObjArray *originalTracks = tracker->GetOriginal();
254 : //
255 : //
256 8 : TTreeSRedirector &cstream = *(tracker->GetDebugStreamer());
257 8 : Int_t ntracks = event->GetNumberOfTracks();
258 8 : Int_t nitstracks = originalTracks->GetEntriesFast();
259 8 : originalTracks->Expand(ntracks);
260 8 : trackHypothesys->Expand(ntracks);
261 8 : bestHypothesys->Expand(ntracks);
262 : //
263 320 : AliHelix * helixes = new AliHelix[ntracks+2];
264 8 : TObjArray trackarray(ntracks+2); //array with tracks - with vertex constrain
265 8 : TObjArray trackarrayc(ntracks+2); //array of "best tracks" - without vertex constrain
266 8 : TObjArray trackarrayl(ntracks+2); //array of "longest tracks" - without vertex constrain
267 : //Change from Bool_t to Int_t for optimization
268 : // Int_t forbN=0;
269 : // Int_t * forbidden = new Int_t [ntracks+2];
270 8 : Bool_t * forbidden = new Bool_t [ntracks+2];
271 16 : Int_t *itsmap = new Int_t [ntracks+2];
272 16 : Float_t *dist = new Float_t[ntracks+2];
273 16 : Float_t *normdist0 = new Float_t[ntracks+2];
274 16 : Float_t *normdist1 = new Float_t[ntracks+2];
275 16 : Float_t *normdist = new Float_t[ntracks+2];
276 16 : Float_t *norm = new Float_t[ntracks+2];
277 16 : Float_t *maxr = new Float_t[ntracks+2];
278 16 : Float_t *minr = new Float_t[ntracks+2];
279 16 : Float_t *minPointAngle= new Float_t[ntracks+2];
280 : //
281 16 : AliV0 *pvertex = new AliV0;
282 16 : AliITStrackMI * dummy= new AliITStrackMI;
283 8 : dummy->SetLabel(0);
284 8 : AliITStrackMI trackat0; //temporary track for DCA calculation
285 : //
286 8 : Float_t primvertex[3]={static_cast<Float_t>(tracker->GetX()),static_cast<Float_t>(tracker->GetY()),static_cast<Float_t>(tracker->GetZ())};
287 : //
288 : // make ITS - ESD map
289 : //
290 320 : for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
291 152 : itsmap[itrack] = -1;
292 : // forbidden[itrack] = 0;
293 152 : forbidden[itrack] = kFALSE;
294 152 : maxr[itrack] = kMaxR;
295 152 : minr[itrack] = kMinR;
296 152 : minPointAngle[itrack] = kMinPointAngle;
297 : }
298 236 : for (Int_t itrack=0;itrack<nitstracks;itrack++){
299 220 : AliITStrackMI * original = (AliITStrackMI*)(originalTracks->At(itrack));
300 110 : Int_t esdindex = original->GetESDtrack()->GetID();
301 110 : itsmap[esdindex] = itrack;
302 : }
303 : //
304 : // create ITS tracks from ESD tracks if not done before
305 : //
306 288 : for (Int_t itrack=0;itrack<ntracks;itrack++){
307 136 : if (itsmap[itrack]>=0) continue;
308 78 : AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
309 : //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
310 : //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
311 26 : tpctrack->GetDZ(tracker->GetX(),tracker->GetY(),tracker->GetZ(),tpctrack->GetDP()); //I.B.
312 46 : if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
313 : // tracks which can reach inner part of ITS
314 : // propagate track to outer its volume - with correction for material
315 16 : AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(tpctrack);
316 : }
317 26 : itsmap[itrack] = nitstracks;
318 26 : originalTracks->AddAt(tpctrack,nitstracks);
319 26 : nitstracks++;
320 26 : }
321 : //
322 : // fill temporary arrays
323 : //
324 288 : for (Int_t itrack=0;itrack<ntracks;itrack++){
325 136 : AliESDtrack * esdtrack = event->GetTrack(itrack);
326 136 : Int_t itsindex = itsmap[itrack];
327 272 : AliITStrackMI *original = (AliITStrackMI*)originalTracks->At(itsmap[itrack]);
328 136 : if (!original) continue;
329 : AliITStrackMI *bestConst = 0;
330 : AliITStrackMI *bestLong = 0;
331 : AliITStrackMI *best = 0;
332 : //
333 : //
334 272 : TObjArray * array = (TObjArray*) trackHypothesys->At(itsindex);
335 354 : Int_t hentries = (array==0) ? 0 : array->GetEntriesFast();
336 : // Get best track with vertex constrain
337 898 : for (Int_t ih=0;ih<hentries;ih++){
338 580 : AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
339 470 : if (!trackh->GetConstrain()) continue;
340 186 : if (!bestConst) bestConst = trackh;
341 220 : if (trackh->GetNumberOfClusters()>kMinClFullTrk ){
342 : bestConst = trackh; // full track - with minimal chi2
343 30 : break;
344 : }
345 320 : if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()) continue;
346 : bestConst = trackh;
347 0 : break;
348 : }
349 : // Get best long track without vertex constrain and best track without vertex constrain
350 966 : for (Int_t ih=0;ih<hentries;ih++){
351 648 : AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
352 438 : if (trackh->GetConstrain()) continue;
353 292 : if (!best) best = trackh;
354 292 : if (!bestLong) bestLong = trackh;
355 420 : if (trackh->GetNumberOfClusters()>kMinClFullTrk){
356 : bestLong = trackh; // full track - with minimal chi2
357 30 : break;
358 : }
359 720 : if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()) continue;
360 : bestLong = trackh;
361 0 : }
362 136 : if (!best) {
363 : best = original;
364 : bestLong = original;
365 54 : }
366 : //I.B. trackat0 = *bestLong;
367 272 : new (&trackat0) AliITStrackMI(*bestLong);
368 : // Double_t xx,yy,zz,alpha;
369 : // if (!bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz)) continue;
370 : // alpha = TMath::ATan2(yy,xx);
371 136 : double alpha = bestLong->PhiPos(); // RS use faster method
372 :
373 : // if (!trackat0.Propagate(alpha,0)) continue;
374 : // trackat0.Propagate(alpha,0); //PH The check on the return value is temporarily disabled (bug 45751)
375 136 : if(!trackat0.Propagate(alpha,0) && kCheckPropagate)continue;
376 : // calculate normalized distances to the vertex
377 : //
378 272 : Float_t ptfac = (1.+100.*TMath::Abs(trackat0.GetC()));
379 272 : if ( bestLong->GetNumberOfClusters()>3 ){
380 156 : dist[itsindex] = trackat0.GetY();
381 78 : norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
382 156 : normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
383 156 : normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
384 78 : normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
385 78 : if (!bestConst){
386 6 : if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
387 6 : if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
388 4 : if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
389 : }else{
390 154 : if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
391 152 : if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
392 : }
393 : }
394 : else{
395 58 : if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
396 0 : dist[itsindex] = bestConst->GetD(0);
397 0 : norm[itsindex] = bestConst->GetDnorm(0);
398 0 : normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
399 0 : normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
400 0 : normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
401 0 : }else{
402 116 : dist[itsindex] = trackat0.GetY();
403 58 : norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
404 116 : normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
405 116 : normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
406 58 : normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
407 116 : if (TMath::Abs(trackat0.GetTgl())>kMinTgl0){
408 2 : if (normdist[itsindex]<kMinNormDistForbTgl0){
409 : // forbN=1;
410 : // forbidden[itsindex]+=1<<forbN;
411 0 : forbidden[itsindex]=kTRUE;
412 0 : }
413 2 : if (normdist[itsindex]>kMinNormDistForbTgl0) {
414 2 : minr[itsindex] = TMath::Max(kMinRTgl0,minr[itsindex]);
415 2 : }
416 : }
417 : }
418 : }
419 :
420 :
421 : //
422 : //-----------------------------------------------------------
423 : //Forbid primary track candidates -
424 : //
425 : //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
426 : //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
427 : //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
428 : //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
429 : //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
430 : //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
431 : //-----------------------------------------------------------
432 136 : if (bestConst){
433 152 : if (bestLong->GetNumberOfClusters()<4 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>kMinClForb0){
434 : // forbN=2;
435 : // forbidden[itsindex]+=1<<forbN;
436 0 : forbidden[itsindex]=kTRUE;
437 0 : }
438 216 : if (normdist[itsindex]<kMinNormDistForb1 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5){
439 : // forbN=3;
440 : // forbidden[itsindex]+=1<<forbN;
441 70 : forbidden[itsindex]=kTRUE;
442 70 : }
443 202 : if (normdist[itsindex]<kMinNormDistForb2 && bestConst->GetClIndex(0)>=0 && bestConst->GetClIndex(1)>=0 ){
444 : // forbN=4;
445 : // forbidden[itsindex]+=1<<forbN;
446 32 : forbidden[itsindex]=kTRUE;
447 32 : }
448 112 : if (normdist[itsindex]<kMinNormDistForb3 && bestConst->GetClIndex(0)>=0){
449 : // forbN=5;
450 : // forbidden[itsindex]+=1<<forbN;
451 32 : forbidden[itsindex]=kTRUE;
452 32 : }
453 150 : if (normdist[itsindex]<kMinNormDistForb4 && bestConst->GetNormChi2(0)<2){
454 : // forbN=6;
455 : // forbidden[itsindex]+=1<<forbN;
456 42 : forbidden[itsindex]=kTRUE;
457 42 : }
458 150 : if (normdist[itsindex]<kMinNormDistForb5 && bestConst->GetNormChi2(0)<1){
459 : // forbN=7;
460 : // forbidden[itsindex]+=1<<forbN;
461 8 : forbidden[itsindex]=kTRUE;
462 8 : }
463 76 : if (bestConst->GetNormChi2(0)<2.5) {
464 56 : minPointAngle[itsindex]= kMinPABestConst;
465 56 : maxr[itsindex] = kMaxRBestConst;
466 56 : }
467 : }
468 : //
469 : //forbid daughter kink candidates
470 : //
471 276 : if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
472 : Bool_t isElectron = kTRUE;
473 : Bool_t isProton = kTRUE;
474 136 : Double_t pid[5];
475 136 : esdtrack->GetESDpid(pid);
476 1360 : for (Int_t i=1;i<5;i++){
477 544 : if (pid[0]<pid[i]) isElectron= kFALSE;
478 544 : if (pid[4]<pid[i]) isProton= kFALSE;
479 : }
480 136 : if (isElectron){
481 136 : forbidden[itsindex]=kFALSE;
482 136 : normdist[itsindex]*=-1;
483 136 : }
484 136 : if (isProton){
485 136 : if (normdist[itsindex]>kMinNormDistForbProt) forbidden[itsindex]=kFALSE;
486 136 : normdist[itsindex]*=-1;
487 136 : }
488 :
489 : // We allow all tracks that are not pions
490 136 : if( (pid[1]+pid[2])< kMaxPidProbPionForb ){
491 136 : forbidden[itsindex]=kFALSE;
492 136 : }
493 :
494 : //
495 : // Causality cuts in TPC volume
496 : //
497 384 : if (esdtrack->GetTPCdensity(0,10) >kMinTPCdensity) maxr[itsindex] = TMath::Min(Float_t(kMaxRTPCdensity0),maxr[itsindex]);
498 394 : if (esdtrack->GetTPCdensity(10,30)>kMinTPCdensity) maxr[itsindex] = TMath::Min(Float_t(kMaxRTPCdensity10),maxr[itsindex]);
499 394 : if (esdtrack->GetTPCdensity(20,40)>kMinTPCdensity) maxr[itsindex] = TMath::Min(Float_t(kMaxRTPCdensity20),maxr[itsindex]);
500 402 : if (esdtrack->GetTPCdensity(30,50)>kMinTPCdensity) maxr[itsindex] = TMath::Min(Float_t(kMaxRTPCdensity30),maxr[itsindex]);
501 : //
502 290 : if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=kMinRTPCdensity;
503 : //
504 : //
505 :
506 272 : if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0){
507 0 : cstream<<"Track"<<
508 0 : "Tr0.="<<best<<
509 0 : "Tr1.="<<((bestConst)? bestConst:dummy)<<
510 0 : "Tr2.="<<bestLong<<
511 0 : "Tr3.="<<&trackat0<<
512 0 : "Esd.="<<esdtrack<<
513 0 : "Dist="<<dist[itsindex]<<
514 0 : "ND0="<<normdist0[itsindex]<<
515 0 : "ND1="<<normdist1[itsindex]<<
516 0 : "ND="<<normdist[itsindex]<<
517 0 : "Pz="<<primvertex[2]<<
518 0 : "Forbid="<<forbidden[itsindex]<<
519 : "\n";
520 : //
521 : }
522 136 : trackarray.AddAt(best,itsindex);
523 136 : trackarrayc.AddAt(bestConst,itsindex);
524 136 : trackarrayl.AddAt(bestLong,itsindex);
525 272 : new (&helixes[itsindex]) AliHelix(*best);
526 136 : }
527 : //
528 : //
529 : //
530 : // first iterration of V0 finder
531 : //
532 : // AM Comment out for optimization
533 : // Int_t rejectBase=0;
534 : // Int_t cutN=0;
535 :
536 288 : for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
537 136 : Int_t itrack0 = itsmap[iesd0];
538 : //-AM comment for optimization and store the forbidden value in the debug streamer
539 136 : if (forbidden[itrack0]) continue;
540 272 : AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
541 136 : if (!btrack0) continue;
542 272 : AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
543 : //
544 2680 : for (Int_t iesd1=iesd0+1;iesd1<ntracks;iesd1++){
545 1204 : Int_t itrack1 = itsmap[iesd1];
546 : //-AM comment for optimization and store the forbidden value in the debug streamer
547 1204 : if (forbidden[itrack1]) continue;
548 :
549 2408 : AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
550 1204 : if (!btrack1) continue;
551 :
552 5328 : if ( (btrack0->GetSign()==btrack1->GetSign()) && !AliITSReconstructor::GetRecoParam()->GetStoreLikeSignV0s()) continue;
553 :
554 : Bool_t isGold = kFALSE;
555 1896 : if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
556 : isGold = kTRUE;
557 28 : }
558 : // rejectBase=0;
559 1264 : AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
560 632 : AliHelix &h1 = helixes[itrack0];
561 632 : AliHelix &h2 = helixes[itrack1];
562 : //
563 : // find linear distance
564 632 : Double_t rmin =0;
565 : //
566 : //
567 : //
568 632 : Double_t phase[2][2]={{0.,0.},{0.,0.}};
569 632 : Double_t radius[2]={0.,0.};
570 632 : Int_t points = h1.GetRPHIintersections(h2, phase, radius);
571 632 : if (points==0) {
572 : // cutN=1;
573 : // rejectBase+=1<<cutN;
574 82 : continue;
575 : }
576 550 : Double_t delta[2]={1000000,1000000};
577 550 : rmin = radius[0];
578 550 : h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
579 550 : if (points==2){
580 720 : if (radius[1]<rmin) rmin = radius[1];
581 486 : h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
582 : }
583 550 : rmin = TMath::Sqrt(rmin);
584 : Double_t distance = 0;
585 : Double_t radiusC = 0;
586 : Int_t iphase = 0;
587 1036 : if (points==1 || delta[0]<delta[1]){
588 294 : distance = TMath::Sqrt(delta[0]);
589 294 : radiusC = TMath::Sqrt(radius[0]);
590 294 : }else{
591 256 : distance = TMath::Sqrt(delta[1]);
592 256 : radiusC = TMath::Sqrt(radius[1]);
593 : iphase=1;
594 : }
595 550 : if (radiusC<TMath::Max(minr[itrack0],minr[itrack1])){
596 : // cutN=2;
597 : //rejectBase+=1<<cutN;
598 204 : continue;
599 : }
600 346 : if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1])){
601 : // cutN=3;
602 : // rejectBase+=1<<cutN;
603 106 : continue;
604 : }
605 240 : Float_t maxDist = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
606 240 : if (distance>maxDist){
607 : // cutN=4;
608 : // rejectBase+=1<<cutN;
609 162 : continue;
610 : }
611 156 : Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
612 78 : if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) {
613 : // cutN=5;
614 : // rejectBase+=1<<cutN;
615 64 : continue;
616 : }
617 : //
618 : //
619 : // Double_t distance = TestV0(h1,h2,pvertex,rmin);
620 : //
621 : // if (distance>maxDist) continue;
622 : // if (pvertex->GetRr()<kMinR) continue;
623 : // if (pvertex->GetRr()>kMaxR) continue;
624 : AliITStrackMI * track0=btrack0;
625 : AliITStrackMI * track1=btrack1;
626 : // if (pvertex->GetRr()<3.5){
627 14 : if (radiusC<3.5){
628 : //use longest tracks inside the pipe
629 0 : track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
630 0 : track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
631 0 : }
632 : //
633 : //
634 :
635 :
636 :
637 14 : pvertex->SetParamN(*track0);
638 14 : pvertex->SetParamP(*track1);
639 14 : pvertex->Update(primvertex);
640 14 : pvertex->SetClusters(track0->ClIndex(),track1->ClIndex()); // register clusters
641 :
642 : // Define gamma, K0, lambda and lambda_bar from the decay particles and calculate the chi2
643 14 : AliKFVertex vertexKF;
644 :
645 28 : vertexKF.X() = tracker->GetX();
646 28 : vertexKF.Y() = tracker->GetY();
647 28 : vertexKF.Z() = tracker->GetZ();
648 28 : vertexKF.Covariance(0,0) = tracker->GetSigmaX()*tracker->GetSigmaX();
649 28 : vertexKF.Covariance(1,1) = tracker->GetSigmaY()*tracker->GetSigmaY();
650 28 : vertexKF.Covariance(2,2) = tracker->GetSigmaZ()*tracker->GetSigmaZ();
651 :
652 14 : AliKFParticle elecKF( *(pvertex->GetParamN()) ,11);
653 14 : AliKFParticle posiKF( *(pvertex->GetParamP()) ,-11);
654 14 : AliKFParticle pipKF( *(pvertex->GetParamN()) , 211);
655 14 : AliKFParticle pinKF( *(pvertex->GetParamP()) ,-211);
656 14 : AliKFParticle protKF( *(pvertex->GetParamP()) ,2212);
657 14 : AliKFParticle aproKF ( *(pvertex->GetParamN()) ,-2212);
658 :
659 :
660 : // Gamma
661 14 : AliKFParticle gamKF(elecKF,posiKF);
662 14 : gamKF.SetProductionVertex(vertexKF);
663 :
664 14 : Float_t gamKFchi2 = 1000;
665 28 : if ( gamKF.GetNDF()!=0 ){
666 42 : gamKFchi2 = gamKF.GetChi2()/gamKF.GetNDF();
667 14 : }
668 :
669 : Float_t massG=0.;
670 : Float_t sigmaMG=0.001;
671 14 : gamKF.SetMassConstraint(massG,sigmaMG);
672 :
673 14 : Float_t gamKFchi2C = 1000;
674 28 : if ( gamKF.GetNDF()!=0 ){
675 42 : gamKFchi2C = gamKF.GetChi2()/gamKF.GetNDF();
676 14 : }
677 :
678 : //K0 short
679 14 : AliKFParticle k0KF(pipKF,pinKF);
680 14 : k0KF.SetProductionVertex(vertexKF);
681 :
682 14 : Float_t k0KFchi2 = 1000;
683 28 : if ( k0KF.GetNDF()!=0 ){
684 42 : k0KFchi2 = k0KF.GetChi2()/k0KF.GetNDF();
685 14 : }
686 :
687 : //Lambda
688 14 : AliKFParticle lambdaKF(protKF,pinKF);
689 14 : lambdaKF.SetProductionVertex(vertexKF);
690 :
691 14 : Float_t lambdaKFchi2 = 1000;
692 28 : if ( lambdaKF.GetNDF()!=0 ){
693 42 : lambdaKFchi2 = lambdaKF.GetChi2()/lambdaKF.GetNDF();
694 14 : }
695 :
696 : //Lambda_bar
697 14 : AliKFParticle alambKF(aproKF,pipKF);
698 14 : alambKF.SetProductionVertex(vertexKF);
699 :
700 14 : Float_t alambKFchi2 = 1000;
701 28 : if ( alambKF.GetNDF()!=0 ){
702 42 : alambKFchi2 = alambKF.GetChi2()/alambKF.GetNDF();
703 14 : }
704 :
705 :
706 :
707 :
708 :
709 14 : if (pvertex->GetRr()<kMinR){
710 : // cutN=6;
711 : // rejectBase+=1<<cutN;
712 0 : continue;
713 : }
714 14 : if (pvertex->GetRr()>kMaxR){
715 : // cutN=7;
716 : // rejectBase+=1<<cutN;
717 0 : continue;
718 : }
719 14 : if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle){
720 : // cutN=8;
721 : // rejectBase+=1<<cutN;
722 0 : continue;
723 : }
724 : //Bo: if (pvertex->GetDist2()>maxDist) continue;
725 :
726 14 : if (pvertex->GetDcaV0Daughters()>maxDist){
727 : // cutN=9;
728 : // rejectBase+=1<<cutN;
729 0 : continue;
730 : }
731 : //Bo: pvertex->SetLab(0,track0->GetLabel());
732 : //Bo: pvertex->SetLab(1,track1->GetLabel());
733 28 : pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
734 28 : pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
735 : //
736 14 : AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
737 14 : AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
738 :
739 : //
740 : //
741 28 : TObjArray * array0b = (TObjArray*)bestHypothesys->At(itrack0);
742 14 : if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<kMinTgl1) {
743 0 : tracker->SetCurrentEsdTrack(itrack0);
744 0 : tracker->FollowProlongationTree((AliITStrackMI*)originalTracks->At(itrack0),itrack0, kFALSE);
745 : }
746 28 : TObjArray * array1b = (TObjArray*)bestHypothesys->At(itrack1);
747 14 : if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<kMinTgl1) {
748 0 : tracker->SetCurrentEsdTrack(itrack1);
749 0 : tracker->FollowProlongationTree((AliITStrackMI*)originalTracks->At(itrack1),itrack1, kFALSE);
750 : }
751 : //
752 28 : AliITStrackMI * track0b = (AliITStrackMI*)originalTracks->At(itrack0);
753 28 : AliITStrackMI * track1b = (AliITStrackMI*)originalTracks->At(itrack1);
754 28 : AliITStrackMI * track0l = (AliITStrackMI*)originalTracks->At(itrack0);
755 28 : AliITStrackMI * track1l = (AliITStrackMI*)originalTracks->At(itrack1);
756 :
757 : Float_t minchi2before0=kMinchi2before0;
758 : Float_t minchi2before1=kMinchi2before1;
759 : Float_t minchi2after0 =kMinchi2after0;
760 : Float_t minchi2after1 =kMinchi2after1;
761 14 : Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
762 14 : Int_t maxLayer = tracker->GetNearestLayer(xrp); //I.B.
763 :
764 176 : if (array0b) for (Int_t i=0;i<5;i++){
765 : // best track after vertex
766 132 : AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
767 102 : if (!btrack) continue;
768 98 : if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
769 : // if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
770 60 : if (btrack->GetX()<pvertex->GetRr()-2.) {
771 58 : if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
772 : Float_t sumchi2= 0;
773 : Float_t sumn = 0;
774 4 : if (maxLayer<3){ // take prim vertex as additional measurement
775 0 : if (normdist[itrack0]>0 && htrackc0){
776 0 : sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
777 0 : }else{
778 0 : sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
779 : }
780 0 : sumn += 3-maxLayer;
781 0 : }
782 36 : for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
783 14 : sumn+=1.;
784 14 : if (btrack->GetClIndex(ilayer)<0){
785 0 : sumchi2+=kAddchi2NegCl0;
786 0 : continue;
787 : }else{
788 14 : Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
789 136 : for (Int_t itrack=0;itrack<4;itrack++){
790 50 : AliITStrackerMI::AliITSlayer &layer=tracker->GetLayer(ilayer);
791 60 : if (layer.GetClusterTracks(itrack,c)>=0 && layer.GetClusterTracks(itrack,c)!=itrack0){
792 2 : sumchi2+=kAddchi2SharedCl; //shared cluster
793 2 : break;
794 : }
795 48 : }
796 14 : sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
797 14 : sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
798 : }
799 14 : }
800 4 : sumchi2/=sumn;
801 8 : if (sumchi2<minchi2before0) minchi2before0=sumchi2;
802 4 : }
803 26 : continue; //safety space - Geo manager will give exact layer
804 : }
805 : track0b = btrack;
806 4 : minchi2after0 = btrack->GetNormChi2(i);
807 4 : break;
808 14 : }
809 170 : if (array1b) for (Int_t i=0;i<5;i++){
810 : // best track after vertex
811 132 : AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
812 112 : if (!btrack) continue;
813 70 : if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
814 : // if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
815 40 : if (btrack->GetX()<pvertex->GetRr()-2){
816 26 : if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
817 : Float_t sumchi2= 0;
818 : Float_t sumn = 0;
819 0 : if (maxLayer<3){ // take prim vertex as additional measurement
820 0 : if (normdist[itrack1]>0 && htrackc1){
821 0 : sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
822 0 : }else{
823 0 : sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
824 : }
825 0 : sumn += 3-maxLayer;
826 0 : }
827 0 : for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
828 0 : sumn+=1.;
829 0 : if (btrack->GetClIndex(ilayer)<0){
830 0 : sumchi2+=kAddchi2NegCl1;
831 0 : continue;
832 : }else{
833 0 : Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
834 0 : for (Int_t itrack=0;itrack<4;itrack++){
835 0 : AliITStrackerMI::AliITSlayer &layer=tracker->GetLayer(ilayer);
836 0 : if (layer.GetClusterTracks(itrack,c)>=0 && layer.GetClusterTracks(itrack,c)!=itrack1){
837 0 : sumchi2+=kAddchi2SharedCl; //shared cluster
838 0 : break;
839 : }
840 0 : }
841 0 : sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
842 0 : sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
843 : }
844 0 : }
845 0 : sumchi2/=sumn;
846 0 : if (sumchi2<minchi2before1) minchi2before1=sumchi2;
847 0 : }
848 14 : continue; //safety space - Geo manager will give exact layer
849 : }
850 : track1b = btrack;
851 6 : minchi2after1 = btrack->GetNormChi2(i);
852 6 : break;
853 14 : }
854 : //
855 : // position resolution - used for DCA cut
856 42 : Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
857 70 : (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
858 56 : (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
859 14 : sigmad =TMath::Sqrt(sigmad)+0.04;
860 14 : if (pvertex->GetRr()>50){
861 6 : Double_t cov0[15],cov1[15];
862 6 : track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
863 6 : track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
864 18 : sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
865 12 : (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
866 6 : (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
867 6 : sigmad =TMath::Sqrt(sigmad)+0.05;
868 6 : }
869 : //
870 14 : AliV0 vertex2;
871 14 : vertex2.SetParamN(*track0b);
872 14 : vertex2.SetParamP(*track1b);
873 14 : vertex2.Update(primvertex);
874 : //Bo: if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
875 26 : if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
876 10 : pvertex->SetParamN(*track0b);
877 10 : pvertex->SetParamP(*track1b);
878 10 : pvertex->Update(primvertex);
879 10 : pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex()); // register clusters
880 20 : pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
881 20 : pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
882 10 : }
883 14 : pvertex->SetDistSigma(sigmad);
884 : //Bo: pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
885 14 : pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
886 : //
887 : // define likelihhod and causalities
888 : //
889 : Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
890 14 : if (maxLayer<1){
891 0 : Float_t fnorm0 = normdist[itrack0];
892 0 : if (fnorm0<0) fnorm0*=-3;
893 0 : Float_t fnorm1 = normdist[itrack1];
894 0 : if (fnorm1<0) fnorm1*=-3;
895 0 : if ((pvertex->GetAnglep()[2]>0.1) || ( (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 ) || (pvertex->GetRr()<3)){
896 0 : pb0 = TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
897 0 : pb1 = TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
898 0 : }
899 0 : pvertex->SetChi2Before(normdist[itrack0]);
900 0 : pvertex->SetChi2After(normdist[itrack1]);
901 0 : pvertex->SetNAfter(0);
902 0 : pvertex->SetNBefore(0);
903 0 : }else{
904 14 : pvertex->SetChi2Before(minchi2before0);
905 14 : pvertex->SetChi2After(minchi2before1);
906 34 : if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
907 4 : pb0 = TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
908 4 : pb1 = TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
909 4 : }
910 14 : pvertex->SetNAfter(maxLayer);
911 14 : pvertex->SetNBefore(maxLayer);
912 : }
913 14 : if (pvertex->GetRr()<90){
914 20 : pa0 *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
915 20 : pa1 *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
916 10 : }
917 14 : if (pvertex->GetRr()<20){
918 2 : pa0 *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
919 2 : pa1 *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
920 2 : }
921 : //
922 14 : pvertex->SetCausality(pb0,pb1,pa0,pa1);
923 : //
924 : // Likelihood calculations - apply cuts
925 : //
926 :
927 : // AM - Comment out for optimization and store the v0OK value
928 : // Int_t v0OK = 0;
929 : // Int_t cutOKN=0;
930 14 : Bool_t v0OK = kTRUE;
931 :
932 :
933 42 : Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
934 42 : p12 += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
935 14 : p12 = TMath::Sqrt(p12); // "mean" momenta
936 :
937 14 : Float_t sigmap0 = kSigp0Par0+kSigp0Par1/(kSigp0Par2+pvertex->GetRr());
938 14 : Float_t sigmap = kSigpPar0*sigmap0*(kSigpPar1+kSigpPar2*p12); // "resolution: of point angle - as a function of radius and momenta
939 :
940 :
941 14 : Float_t causalityA = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
942 28 : Float_t causalityB = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Double_t(0.7))*
943 14 : TMath::Min(pvertex->GetCausalityP()[3],Double_t(0.7)));
944 : //
945 : //Bo: Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
946 14 : Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
947 14 : Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<kMaxDcaLh0)*(lDistNorm<5);
948 :
949 42 : Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
950 28 : 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
951 28 : 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
952 14 : 0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
953 : //
954 14 : if (causalityA<kCausality0Cut){
955 : // cutOKN=1;
956 : // v0OK += 1<<cutOKN;
957 2 : v0OK = kFALSE;
958 2 : }
959 14 : if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut){
960 : // cutOKN=2;
961 : // v0OK += 1<<cutOKN;
962 2 : v0OK = kFALSE;
963 2 : }
964 14 : if (likelihood1<kLikelihood1Cut){
965 : // cutOKN=3;
966 : // v0OK += 1<<cutOKN;
967 2 : v0OK = kFALSE;
968 2 : }
969 14 : if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut){
970 : // cutOKN=4;
971 : // v0OK += 1<<cutOKN;
972 2 : v0OK = kFALSE;
973 2 : }
974 : //
975 : //
976 28 : if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0){
977 0 : Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
978 0 : cstream<<"It0"<<
979 0 : "Tr0.="<<track0<< //best without constrain
980 0 : "Tr1.="<<track1<< //best without constrain
981 0 : "posiKF.="<<&posiKF<< //KF from track0
982 0 : "elecKF.="<<&elecKF<< //KF from track1
983 0 : "Tr0B.="<<track0b<< //best without constrain after vertex
984 0 : "Tr1B.="<<track1b<< //best without constrain after vertex
985 0 : "Tr0C.="<<htrackc0<< //best with constrain if exist
986 0 : "Tr1C.="<<htrackc1<< //best with constrain if exist
987 0 : "Tr0L.="<<track0l<< //longest best
988 0 : "Tr1L.="<<track1l<< //longest best
989 0 : "Esd0.="<<track0->GetESDtrack()<< // esd track0 params
990 0 : "Esd1.="<<track1->GetESDtrack()<< // esd track1 params
991 0 : "V0.="<<pvertex<< //vertex properties
992 0 : "V0b.="<<&vertex2<< //vertex properties at "best" track
993 0 : "gamKF.="<<&gamKF<< //KF from pvertex
994 0 : "k0KF.="<<&k0KF<< //KF from pvertex
995 0 : "lambdaKF.="<<&lambdaKF<< //KF from pvertex
996 0 : "alambKF.="<<&lambdaKF<< //KF from pvertex
997 0 : "gamKFchi2="<<gamKFchi2<< //Normalized chi2 from KF assuming gamma momther
998 0 : "gamKFchi2C="<<gamKFchi2C<< //Normalized chi2 from KF assuming gamma mother+mass constraint
999 0 : "k0KFchi2="<<k0KFchi2<< //Normalized chi2 from KF assuming K0 mother
1000 0 : "lambdaKFchi2="<<lambdaKFchi2<< //Normalized chi2 from KF assuming Lambda mother
1001 0 : "alambKFchi2="<<alambKFchi2<< //Normalized chi2 from KF assuming lambda_bar mother
1002 0 : "ND0="<<normdist[itrack0]<< //normalize distance for track0
1003 0 : "ND1="<<normdist[itrack1]<< //normalize distance for track1
1004 0 : "Gold="<<gold<< //
1005 : // "RejectBase="<<rejectBase<< //rejection in First itteration
1006 0 : "OK="<<v0OK<<
1007 0 : "rmin="<<rmin<<
1008 0 : "sigmad="<<sigmad<<
1009 0 : "Forbid0="<<forbidden[itrack0]<<
1010 0 : "Forbid1="<<forbidden[itrack1]<<
1011 : "\n";
1012 0 : }
1013 : // if (rejectBase>0) continue;
1014 : //if (forbidden[itrack0]>0 ||forbidden[itrack1]>0) continue;
1015 : //
1016 14 : pvertex->SetStatus(0);
1017 : // if (rejectBase) {
1018 : // pvertex->SetStatus(-100);
1019 : //}
1020 14 : if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
1021 : //Bo: pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
1022 28 : pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
1023 28 : pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
1024 : // if (v0OK==0){
1025 14 : if (v0OK){
1026 : // AliV0vertex vertexjuri(*track0,*track1);
1027 : // vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
1028 : // event->AddV0(&vertexjuri);
1029 10 : pvertex->SetStatus(100);
1030 10 : }
1031 14 : pvertex->SetOnFlyStatus(kTRUE);
1032 14 : pvertex->ChangeMassHypothesis(kK0Short);
1033 14 : event->AddV0(pvertex);
1034 : }
1035 1210 : }
1036 136 : }
1037 : //
1038 : //
1039 : // delete temporary arrays
1040 : //
1041 16 : delete[] forbidden;
1042 16 : delete[] minPointAngle;
1043 16 : delete[] maxr;
1044 16 : delete[] minr;
1045 16 : delete[] norm;
1046 16 : delete[] normdist;
1047 16 : delete[] normdist1;
1048 16 : delete[] normdist0;
1049 16 : delete[] dist;
1050 16 : delete[] itsmap;
1051 176 : delete[] helixes;
1052 16 : delete pvertex;
1053 16 : delete dummy;
1054 16 : }
1055 : //------------------------------------------------------------------------
1056 : void AliITSV0Finder::RefitV02(const AliESDEvent *event,
1057 : AliITStrackerMI *tracker)
1058 : {
1059 : //
1060 : //try to refit V0s in the third path of the reconstruction
1061 : //
1062 16 : TTreeSRedirector &cstream = *(tracker->GetDebugStreamer());
1063 : //
1064 8 : Int_t nv0s = event->GetNumberOfV0s();
1065 8 : Float_t primvertex[3]={static_cast<Float_t>(tracker->GetX()),static_cast<Float_t>(tracker->GetY()),static_cast<Float_t>(tracker->GetZ())};
1066 8 : AliV0 v0temp;
1067 44 : for (Int_t iv0 = 0; iv0<nv0s;iv0++){
1068 28 : AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
1069 14 : if (!v0mi) continue;
1070 14 : Int_t itrack0 = v0mi->GetIndex(0);
1071 14 : Int_t itrack1 = v0mi->GetIndex(1);
1072 14 : AliESDtrack *esd0 = event->GetTrack(itrack0);
1073 14 : AliESDtrack *esd1 = event->GetTrack(itrack1);
1074 14 : if (!esd0||!esd1) continue;
1075 14 : AliITStrackMI tpc0(*esd0);
1076 14 : AliITStrackMI tpc1(*esd1);
1077 14 : Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
1078 14 : Double_t alpha =TMath::ATan2(y,x); //I.B.
1079 14 : if (v0mi->GetRr()>85){
1080 16 : if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
1081 4 : v0temp.SetParamN(tpc0);
1082 4 : v0temp.SetParamP(tpc1);
1083 4 : v0temp.Update(primvertex);
1084 8 : if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0) cstream<<"Refit"<<
1085 0 : "V0.="<<v0mi<<
1086 0 : "V0refit.="<<&v0temp<<
1087 0 : "Tr0.="<<&tpc0<<
1088 0 : "Tr1.="<<&tpc1<<
1089 : "\n";
1090 : //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
1091 8 : if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
1092 2 : v0mi->SetParamN(tpc0);
1093 2 : v0mi->SetParamP(tpc1);
1094 2 : v0mi->Update(primvertex);
1095 : }
1096 : }
1097 4 : continue;
1098 : }
1099 10 : if (v0mi->GetRr()>35){
1100 6 : AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
1101 6 : AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
1102 24 : if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
1103 6 : v0temp.SetParamN(tpc0);
1104 6 : v0temp.SetParamP(tpc1);
1105 6 : v0temp.Update(primvertex);
1106 12 : if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0) cstream<<"Refit"<<
1107 0 : "V0.="<<v0mi<<
1108 0 : "V0refit.="<<&v0temp<<
1109 0 : "Tr0.="<<&tpc0<<
1110 0 : "Tr1.="<<&tpc1<<
1111 : "\n";
1112 : //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
1113 10 : if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
1114 4 : v0mi->SetParamN(tpc0);
1115 4 : v0mi->SetParamP(tpc1);
1116 4 : v0mi->Update(primvertex);
1117 : }
1118 : }
1119 6 : continue;
1120 : }
1121 4 : AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
1122 4 : AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
1123 : // if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
1124 12 : if (tracker->RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) &&
1125 4 : tracker->RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
1126 4 : v0temp.SetParamN(tpc0);
1127 4 : v0temp.SetParamP(tpc1);
1128 4 : v0temp.Update(primvertex);
1129 8 : if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0) cstream<<"Refit"<<
1130 0 : "V0.="<<v0mi<<
1131 0 : "V0refit.="<<&v0temp<<
1132 0 : "Tr0.="<<&tpc0<<
1133 0 : "Tr1.="<<&tpc1<<
1134 : "\n";
1135 : //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
1136 8 : if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
1137 2 : v0mi->SetParamN(tpc0);
1138 2 : v0mi->SetParamP(tpc1);
1139 2 : v0mi->Update(primvertex);
1140 : }
1141 : }
1142 18 : }
1143 8 : }
1144 : //------------------------------------------------------------------------
1145 :
1146 :
1147 : AliITSV0Finder::~AliITSV0Finder()
1148 0 : {
1149 : //
1150 : //destructor
1151 : //
1152 0 : if (fDebugStreamer) {
1153 : //fDebugStreamer->Close();
1154 0 : delete fDebugStreamer;
1155 : }
1156 :
1157 0 : }
|
