Fixing compilation warnings
[u/mrichter/AliRoot.git] / ITS / AliITSVertexer3D.cxx
1 /**************************************************************************
2  * Copyright(c) 2006-2008, 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 #include <TTree.h>
16 #include "AliESDVertex.h"
17 #include "AliLog.h"
18 #include "AliStrLine.h"
19 #include "AliTracker.h"
20 #include "AliITSDetTypeRec.h"
21 #include "AliITSRecPoint.h"
22 #include "AliITSgeomTGeo.h"
23 #include "AliVertexerTracks.h"
24 #include "AliITSVertexer3D.h"
25 #include "AliITSVertexerZ.h"
26 #include "AliITSSortTrkl.h"
27 /////////////////////////////////////////////////////////////////
28 // this class implements a method to determine
29 // the 3 coordinates of the primary vertex
30 // for p-p collisions 
31 // It can be used successfully with Pb-Pb collisions
32 ////////////////////////////////////////////////////////////////
33
34 ClassImp(AliITSVertexer3D)
35
36 /* $Id$ */
37
38 //______________________________________________________________________
39 AliITSVertexer3D::AliITSVertexer3D():
40   AliITSVertexer(),
41   fLines("AliStrLine",1000),
42   fVert3D(),
43   fCoarseDiffPhiCut(0.),
44   fFineDiffPhiCut(0.),
45   fCutOnPairs(0.),
46   fCoarseMaxRCut(0.),
47   fMaxRCut(0.),
48   fMaxRCut2(0.),
49   fZCutDiamond(0.),
50   fMaxZCut(0.),
51   fDCAcut(0.),
52   fDiffPhiMax(0.),
53   fMeanPSelTrk(0.),
54   fMeanPtSelTrk(0.),
55   fUsedCluster(kMaxCluPerMod*kNSPDMod),
56   fZHisto(0),
57   fDCAforPileup(0.),
58   fDiffPhiforPileup(0.),
59   fBinSizeR(0.),
60   fBinSizeZ(0.),
61   fPileupAlgo(0)
62 {
63   // Default constructor
64   SetCoarseDiffPhiCut();
65   SetFineDiffPhiCut();
66   SetCutOnPairs();
67   SetCoarseMaxRCut();
68   SetMaxRCut();
69   SetMaxRCutAlgo2();
70   SetZCutDiamond();
71   SetMaxZCut();
72   SetDCACut();
73   SetDiffPhiMax();
74   SetMeanPSelTracks();
75   SetMeanPtSelTracks();
76   SetMinDCAforPileup();
77   SetDeltaPhiforPileup();
78   SetPileupAlgo();
79   SetBinSizeR();
80   SetBinSizeZ();
81   Double_t binsize=0.02; // default 200 micron
82   Int_t nbins=static_cast<Int_t>(1+2*fZCutDiamond/binsize);
83   fZHisto=new TH1F("hz","",nbins,-fZCutDiamond,-fZCutDiamond+binsize*nbins);
84 }
85
86 //______________________________________________________________________
87 AliITSVertexer3D::~AliITSVertexer3D() {
88   // Destructor
89   fLines.Clear("C");
90   if(fZHisto) delete fZHisto;
91 }
92
93 //______________________________________________________________________
94 void AliITSVertexer3D::ResetVert3D(){
95   //
96   ResetVertex();
97   fVert3D.SetXv(0.);
98   fVert3D.SetYv(0.);
99   fVert3D.SetZv(0.);
100   fVert3D.SetDispersion(0.);
101   fVert3D.SetNContributors(0);
102   fUsedCluster.ResetAllBits(0);
103 }
104 //______________________________________________________________________
105 AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree){
106   // Defines the AliESDVertex for the current event
107   ResetVert3D();
108   AliDebug(1,"FindVertexForCurrentEvent - 3D - PROCESSING NEXT EVENT");
109   fLines.Clear("C");
110   fCurrentVertex = NULL;
111
112   Int_t nolines = FindTracklets(itsClusterTree,0);
113   Int_t rc;
114   if(nolines>=2){
115     rc=Prepare3DVertex(0);
116     if(fVert3D.GetNContributors()>0){
117       fLines.Clear("C");
118       nolines = FindTracklets(itsClusterTree,1);
119       if(nolines>=2){
120         rc=Prepare3DVertex(1);
121         if(fPileupAlgo == 2 && rc == 0) FindVertex3DIterative();
122         else if(fPileupAlgo<2 && rc == 0) FindVertex3D(itsClusterTree);
123         if(rc!=0) fVert3D.SetNContributors(0); // exclude this vertex      
124       }
125     }
126   }
127
128   if(!fCurrentVertex){
129     AliITSVertexerZ vertz(GetNominalPos()[0],GetNominalPos()[1]);
130     vertz.SetDetTypeRec(GetDetTypeRec());
131     AliDebug(1,"Call Vertexer Z\n");
132     vertz.SetLowLimit(-fZCutDiamond);
133     vertz.SetHighLimit(fZCutDiamond);
134     AliESDVertex* vtxz = vertz.FindVertexForCurrentEvent(itsClusterTree);
135     if(vtxz){
136       Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],vtxz->GetZv()};
137       Double_t covmatrix[6];
138       vtxz->GetCovMatrix(covmatrix);
139       Double_t chi2=99999.;
140       Int_t    nContr=vtxz->GetNContributors();
141       fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);    
142       fCurrentVertex->SetTitle("vertexer: Z");
143       fCurrentVertex->SetName("SPDVertexZ");
144       delete vtxz;
145     }
146
147   }
148   FindMultiplicity(itsClusterTree);
149   return fCurrentVertex;
150 }  
151
152 //______________________________________________________________________
153 void AliITSVertexer3D::FindVertex3D(TTree *itsClusterTree){
154  
155   Double_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv());
156   if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){
157     Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
158     Double_t covmatrix[6];
159     fVert3D.GetCovMatrix(covmatrix);
160     Double_t chi2=99999.;
161     Int_t    nContr=fVert3D.GetNContributors();
162     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);    
163     fCurrentVertex->SetTitle("vertexer: 3D");
164     fCurrentVertex->SetName("SPDVertex3D");
165     fCurrentVertex->SetDispersion(fVert3D.GetDispersion());
166     fNoVertices=1;
167     
168     switch(fPileupAlgo){
169     case 0: PileupFromZ(); break;
170     case 1: FindOther3DVertices(itsClusterTree); break;
171     default: AliError("Wrong pileup algorithm"); break;
172     }
173     if(fNoVertices==1){
174       fVertArray = new AliESDVertex[1];
175       fVertArray[0]=(*fCurrentVertex);    
176     }
177   }
178 }
179
180 //______________________________________________________________________
181 void AliITSVertexer3D::FindVertex3DIterative(){
182   //
183
184   Int_t nLines=fLines.GetEntriesFast();
185   Int_t maxPoints=nLines*(nLines-1)/2;
186   Double_t* xP=new Double_t[maxPoints];
187   Double_t* yP=new Double_t[maxPoints];
188   Double_t* zP=new Double_t[maxPoints];
189   Int_t* index1=new Int_t[maxPoints];
190   Int_t* index2=new Int_t[maxPoints];
191   Double_t xbeam=fVert3D.GetXv();
192   Double_t ybeam=fVert3D.GetYv();
193
194   Int_t iPoint=0;
195   for(Int_t ilin1=0; ilin1<nLines; ilin1++){
196     AliStrLine *l1 = (AliStrLine*)fLines.At(ilin1);
197     for(Int_t ilin2=ilin1+1; ilin2<nLines; ilin2++){
198       AliStrLine *l2 = (AliStrLine*)fLines.At(ilin2);
199       Double_t dca=l1->GetDCA(l2);
200       if(dca > fDCAcut || dca<0.00001) continue;
201       Double_t point[3];
202       Int_t retc = l1->Cross(l2,point);
203       if(retc<0)continue;
204       Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
205       if(rad>fCoarseMaxRCut)continue;
206       Double_t distFromBeam=TMath::Sqrt((point[0]-xbeam)*(point[0]-xbeam)+(point[1]-ybeam)*(point[1]-ybeam));
207       if(distFromBeam>fMaxRCut2) continue;
208       xP[iPoint]=point[0];
209       yP[iPoint]=point[1];
210       zP[iPoint]=point[2];
211       index1[iPoint]=ilin1;
212       index2[iPoint]=ilin2;
213       iPoint++;
214     }
215   }
216   Int_t npoints=iPoint++;
217   Int_t index=0;
218   Short_t* mask=new Short_t[npoints];
219   for(Int_t ip=0;ip<npoints;ip++) mask[ip]=-1;
220  
221   for(Int_t ip1=0;ip1<npoints;ip1++){
222     if(mask[ip1]==-1) mask[ip1]=index++;
223     for(Int_t ip2=ip1+1; ip2<npoints; ip2++){
224       if(mask[ip2]==mask[ip1] && mask[ip2]!=-1) continue;
225       Double_t dist2=(xP[ip1]-xP[ip2])*(xP[ip1]-xP[ip2]);
226       dist2+=(yP[ip1]-yP[ip2])*(yP[ip1]-yP[ip2]);
227       dist2+=(zP[ip1]-zP[ip2])*(zP[ip1]-zP[ip2]);
228       if(dist2<fCutOnPairs*fCutOnPairs){ 
229         if(mask[ip2]==-1) mask[ip2]=mask[ip1];
230         else{
231           for(Int_t ip=0; ip<npoints;ip++){
232             if(mask[ip]==mask[ip2]) mask[ip]=mask[ip1];
233           }
234         }
235       }
236     }
237   }
238
239
240   // Count multiplicity of trackelts in clusters
241   UInt_t* isIndUsed=new UInt_t[index+1];
242   for(Int_t ind=0;ind<index+1;ind++) isIndUsed[ind]=0;
243   for(Int_t ip=0; ip<npoints;ip++){
244     Int_t ind=mask[ip];
245     isIndUsed[ind]++;
246   }
247
248   // Count clusters/vertices and sort according to multiplicity
249   Int_t nClusters=0;
250   Int_t* sortedIndex=new Int_t[index+1];
251   for(Int_t ind1=0;ind1<index+1;ind1++){
252     if(isIndUsed[ind1]<=1) isIndUsed[ind1]=0;
253     else nClusters++;
254     UInt_t cap=9999999;
255     if(ind1>0) cap=isIndUsed[sortedIndex[ind1-1]];
256     UInt_t bigger=0;
257     Int_t biggerindex=-1;
258     for(Int_t ind2=0;ind2<index+1;ind2++){
259       Bool_t use=kTRUE;
260       for(Int_t ind3=0; ind3<ind1; ind3++)
261         if(ind2==sortedIndex[ind3]) use=kFALSE;
262       if(use && isIndUsed[ind2]>bigger && isIndUsed[ind2]<=cap){
263         bigger=isIndUsed[ind2];
264         biggerindex=ind2;
265       }
266     }
267     sortedIndex[ind1]=biggerindex;    
268   }
269   AliDebug(3,Form("Number of clusters before merging = %d\n",nClusters));
270
271   // Assign lines to clusters/vertices and merge clusters which share 1 line
272   Int_t nClustersAfterMerge=nClusters;
273   Int_t* belongsTo=new Int_t[nLines];
274   for(Int_t ilin=0; ilin<nLines; ilin++) belongsTo[ilin]=-1;
275   for(Int_t iclu=0;iclu<nClusters;iclu++){
276     Int_t actualCluIndex=iclu;
277     for(Int_t ip=0; ip<npoints;ip++){
278       if(mask[ip]==sortedIndex[iclu]){
279         Int_t ind1=index1[ip];
280         if(belongsTo[ind1]==-1) belongsTo[ind1]=actualCluIndex;
281         else if(belongsTo[ind1]<actualCluIndex){
282           Int_t newCluIndex=belongsTo[ind1];
283           for(Int_t ilin=0; ilin<nLines; ilin++){
284             if(belongsTo[ilin]==actualCluIndex) belongsTo[ilin]=newCluIndex;
285           }
286           AliDebug(10,Form("Merged cluster %d with %d\n",actualCluIndex,newCluIndex));
287           actualCluIndex=newCluIndex;
288           nClustersAfterMerge--;
289         }
290         Int_t ind2=index2[ip];      
291         if(belongsTo[ind2]==-1) belongsTo[ind2]=actualCluIndex;
292         else if(belongsTo[ind2]<actualCluIndex){
293           Int_t newCluIndex=belongsTo[ind2];
294           for(Int_t ilin=0; ilin<nLines; ilin++){
295             if(belongsTo[ilin]==actualCluIndex) belongsTo[ilin]=newCluIndex;
296           }
297           AliDebug(10,Form("Merged cluster %d with %d\n",actualCluIndex,newCluIndex));
298           actualCluIndex=newCluIndex;
299           nClustersAfterMerge--;
300         }
301       }
302     }
303   }
304   AliDebug(3,Form("Number of clusters after merging = %d\n",nClustersAfterMerge));
305   
306   // Count lines associated to each cluster/vertex
307   UInt_t *cluSize=new UInt_t[nClusters];
308   for(Int_t iclu=0;iclu<nClusters;iclu++){ 
309     cluSize[iclu]=0;
310     for(Int_t ilin=0; ilin<nLines; ilin++){
311       if(belongsTo[ilin]==iclu) cluSize[iclu]++;
312     }
313   }
314
315   // Count good vertices (>1 associated tracklet)
316   UInt_t nGoodVert=0;
317   for(Int_t iclu=0;iclu<nClusters;iclu++){ 
318     AliDebug(3,Form("Vertex %d Size=%d\n",iclu,cluSize[iclu]));
319     if(cluSize[iclu]>1) nGoodVert++;
320   }
321     
322   AliDebug(1,Form("Number of good vertices = %d\n",nGoodVert));
323   // Calculate vertex coordinates for each cluster
324   if(nGoodVert>0){
325     fVertArray = new AliESDVertex[nGoodVert];
326     Int_t iVert=0;
327     for(Int_t iclu=0;iclu<nClusters;iclu++){
328       Int_t size=cluSize[iclu];
329       if(size>1){
330         AliStrLine **arrlin = new AliStrLine*[size];
331         Int_t nFilled=0;
332         for(Int_t ilin=0; ilin<nLines; ilin++){
333           if(belongsTo[ilin]==iclu){
334             arrlin[nFilled++] = dynamic_cast<AliStrLine*>(fLines[ilin]);
335           }
336         }      
337         AliDebug(3,Form("Vertex %d  N associated tracklets = %d out of %d\n",iVert,size,nFilled));
338
339         fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin,nFilled);
340         Double_t peak[3];
341         fVertArray[iVert].GetXYZ(peak);
342         AliStrLine **arrlin2 = new AliStrLine*[size];
343         Int_t nFilled2=0;       
344         for(Int_t i=0; i<nFilled;i++){
345           AliStrLine *l1 = arrlin[i];     
346           if(l1->GetDistFromPoint(peak)< fDCAcut)
347             arrlin2[nFilled2++] = dynamic_cast<AliStrLine*>(l1);
348         }
349         if(nFilled2>1){
350           AliDebug(3,Form("Vertex %d  recalculated with %d tracklets\n",iVert,nFilled2));
351           fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin2,nFilled2);
352         }
353         delete [] arrlin;
354         delete [] arrlin2;
355         ++iVert;
356       }
357     }
358     
359     if(nGoodVert > 1){
360       fIsPileup = kTRUE;
361       fNTrpuv = fVertArray[1].GetNContributors();
362       fZpuv = fVertArray[1].GetZv();
363     }
364     
365     Double_t vRadius=TMath::Sqrt(fVertArray[0].GetXv()*fVertArray[0].GetXv()+fVertArray[0].GetYv()*fVertArray[0].GetYv());
366     if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){
367       Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()};
368       Double_t covmatrix[6];
369       fVertArray[0].GetCovMatrix(covmatrix);
370       Double_t chi2=99999.;
371       Int_t    nContr=fVertArray[0].GetNContributors();
372       fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);    
373       fCurrentVertex->SetTitle("vertexer: 3D");
374       fCurrentVertex->SetName("SPDVertex3D");
375       fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion());  
376     }
377   }
378
379   delete [] index1;
380   delete [] index2;
381   delete [] mask;
382   delete [] isIndUsed;
383   delete [] sortedIndex;
384   delete [] belongsTo;
385   delete [] cluSize;
386   delete [] xP;
387   delete [] yP;
388   delete [] zP;
389 }
390 //______________________________________________________________________
391 void AliITSVertexer3D::FindVertex3DIterativeMM(){
392   // Defines the AliESDVertex for the current event
393   Int_t numsor=fLines.GetEntriesFast()*(fLines.GetEntriesFast()-1)/2;
394   //cout<<"AliITSVertexer3D::FindVertexForCurentEvent: Number of tracklets selected for vertexing "<<fLines.GetEntriesFast()<<"; Number of pairs: "<<numsor<<endl;
395   AliITSSortTrkl srt(fLines,numsor,fCutOnPairs,fCoarseMaxRCut); 
396   srt.FindClusters();   
397   AliInfo(Form("Number of vertices: %d",srt.GetNumberOfClusters()));
398       
399   fNoVertices = srt.GetNumberOfClusters();
400   //printf("fNoVertices = %d \n",fNoVertices);
401   if(fNoVertices>0){
402     fVertArray = new AliESDVertex[fNoVertices];
403     for(Int_t kk=0; kk<srt.GetNumberOfClusters(); kk++){
404       Int_t size = 0;
405       Int_t *labels = srt.GetTrackletsLab(kk,size);
406       /*
407         Int_t *pairs = srt.GetClusters(kk);
408         Int_t nopai = srt.GetSizeOfCluster(kk);
409         cout<<"***** Vertex number "<<kk<<".  Pairs: \n";
410         for(Int_t jj=0;jj<nopai;jj++){
411         cout<<pairs[jj]<<" - ";
412         if(jj>0 & jj%8==0)cout<<endl;
413         }
414         cout<<endl;
415         cout<<"***** Vertex number "<<kk<<".  Labels: \n";
416       */
417       AliStrLine **tclo = new AliStrLine* [size];
418       for(Int_t jj=0;jj<size;jj++){
419         //          cout<<labels[jj]<<" - ";
420         //          if(jj>0 & jj%8==0)cout<<endl;
421         tclo[jj] = dynamic_cast<AliStrLine*>(fLines[labels[jj]]);
422       }
423       //          cout<<endl;
424       delete []labels;
425       fVertArray[kk]=AliVertexerTracks::TrackletVertexFinder(tclo,size);
426       delete [] tclo;
427       //          fVertArray[kk].PrintStatus();
428       if(kk == 1){
429         // at least one second vertex is present
430         fIsPileup = kTRUE;
431         fNTrpuv = fVertArray[kk].GetNContributors();
432         fZpuv = fVertArray[kk].GetZv();
433       }
434     }
435     Double_t vRadius=TMath::Sqrt(fVertArray[0].GetXv()*fVertArray[0].GetXv()+fVertArray[0].GetYv()*fVertArray[0].GetYv());
436     if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){
437       Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()};
438       Double_t covmatrix[6];
439       fVertArray[0].GetCovMatrix(covmatrix);
440       Double_t chi2=99999.;
441       Int_t    nContr=fVertArray[0].GetNContributors();
442       fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);    
443       fCurrentVertex->SetTitle("vertexer: 3D");
444       fCurrentVertex->SetName("SPDVertex3D");
445       fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion());  
446     }
447   }
448
449 }  
450
451 //______________________________________________________________________
452 Bool_t AliITSVertexer3D::DistBetweenVertices(AliESDVertex &a, AliESDVertex &b, Double_t test, Double_t &dist){
453   // method to compare the distance between vertices a and b with "test"
454   //it returns kTRUE is the distance is less or equal to test
455   dist = (a.GetX()-b.GetX()) * (a.GetX()-b.GetX());
456   dist +=  (a.GetY()-b.GetY()) * (a.GetY()-b.GetY());
457   dist +=  (a.GetZ()-b.GetZ()) * (a.GetZ()-b.GetZ());
458   dist = TMath::Sqrt(dist);
459   if(dist <= test)return kTRUE;
460   return kFALSE;
461 }
462
463
464 //______________________________________________________________________
465 Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){
466   // All the possible combinations between recpoints on layer 1and 2 are
467   // considered. Straight lines (=tracklets)are formed. 
468   // The tracklets are processed in Prepare3DVertex
469
470   if(!GetDetTypeRec())AliFatal("DetTypeRec pointer has not been set");
471
472   TTree *tR = itsClusterTree;
473   fDetTypeRec->ResetRecPoints();
474   fDetTypeRec->SetTreeAddressR(tR);
475   TClonesArray *itsRec  = 0;
476   if(optCuts==0) fZHisto->Reset();
477  // gc1 are local and global coordinates for layer 1
478   Float_t gc1f[3]={0.,0.,0.};
479   Double_t gc1[3]={0.,0.,0.};
480   // gc2 are local and global coordinates for layer 2
481   Float_t gc2f[3]={0.,0.,0.};
482   Double_t gc2[3]={0.,0.,0.};
483
484   itsRec = fDetTypeRec->RecPoints();
485   TBranch *branch = NULL;
486   branch = tR->GetBranch("ITSRecPoints");
487   if(!branch){
488     AliWarning("Null pointer for RecPoints branch, vertex not calculated");
489     return -1;
490   }
491
492   // Set values for cuts
493   Double_t xbeam=GetNominalPos()[0]; 
494   Double_t ybeam=GetNominalPos()[1];
495   Double_t zvert=0.;
496   Double_t deltaPhi=fCoarseDiffPhiCut;
497   Double_t deltaR=fCoarseMaxRCut;
498   Double_t dZmax=fZCutDiamond;
499   if(optCuts==1){
500     xbeam=fVert3D.GetXv();
501     ybeam=fVert3D.GetYv();
502     zvert=fVert3D.GetZv();
503     deltaPhi = fDiffPhiMax; 
504     deltaR=fMaxRCut;
505     dZmax=fMaxZCut;
506     if(fPileupAlgo == 2){
507       dZmax=fZCutDiamond;
508       deltaR=fMaxRCut2;
509     }
510   } else if(optCuts==2){
511     xbeam=fVert3D.GetXv();
512     ybeam=fVert3D.GetYv();
513     deltaPhi = fDiffPhiforPileup;
514     deltaR=fMaxRCut;
515   }
516
517   Int_t nrpL1 = 0;    // number of rec points on layer 1
518   Int_t nrpL2 = 0;    // number of rec points on layer 2
519
520   // By default irstL1=0 and lastL1=79
521   Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
522   Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;
523   for(Int_t module= firstL1; module<=lastL1;module++){  // count number of recopints on layer 1
524     branch->GetEvent(module);
525     nrpL1+= itsRec->GetEntries();
526     fDetTypeRec->ResetRecPoints();
527   }
528   //By default firstL2=80 and lastL2=239
529   Int_t firstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
530   Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1;
531   for(Int_t module= firstL2; module<=lastL2;module++){  // count number of recopints on layer 2
532     branch->GetEvent(module);
533     nrpL2+= itsRec->GetEntries();
534     fDetTypeRec->ResetRecPoints();
535   }
536   if(nrpL1 == 0 || nrpL2 == 0){
537     AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",nrpL1,nrpL2));
538     return -1;
539   }
540   AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2));
541
542   Double_t a[3]={xbeam,ybeam,0.}; 
543   Double_t b[3]={xbeam,ybeam,10.};
544   AliStrLine zeta(a,b,kTRUE);
545   static Double_t bField=TMath::Abs(AliTracker::GetBz()/10.); //T
546   SetMeanPPtSelTracks(bField);
547
548   Int_t nolines = 0;
549   // Loop on modules of layer 1
550   for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){   // Loop on modules of layer 1
551     if(!fUseModule[modul1]) continue;
552     UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
553     branch->GetEvent(modul1);
554     Int_t nrecp1 = itsRec->GetEntries();
555     static TClonesArray prpl1("AliITSRecPoint",nrecp1);
556     prpl1.SetOwner();
557     for(Int_t j=0;j<nrecp1;j++){
558       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
559       new(prpl1[j])AliITSRecPoint(*recp);
560     }
561     fDetTypeRec->ResetRecPoints();
562     for(Int_t j=0;j<nrecp1;j++){
563       if(j>kMaxCluPerMod) continue;
564       UShort_t idClu1=modul1*kMaxCluPerMod+j;
565       if(fUsedCluster.TestBitNumber(idClu1)) continue;
566       AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1.At(j);
567       recp1->GetGlobalXYZ(gc1f);
568       for(Int_t ico=0;ico<3;ico++)gc1[ico]=gc1f[ico];
569
570       Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam);
571       if(phi1<0)phi1=2*TMath::Pi()+phi1;
572       for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
573         for(Int_t k=0;k<4;k++){
574           Int_t ladmod=fLadders[ladder-1]+ladl2;
575           if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
576           Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
577           if(!fUseModule[modul2]) continue;
578           branch->GetEvent(modul2);
579           Int_t nrecp2 = itsRec->GetEntries();
580           for(Int_t j2=0;j2<nrecp2;j2++){
581             if(j2>kMaxCluPerMod) continue;
582             UShort_t idClu2=modul2*kMaxCluPerMod+j2;
583             if(fUsedCluster.TestBitNumber(idClu2)) continue;
584             AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2);
585             recp2->GetGlobalXYZ(gc2f);
586             for(Int_t ico=0;ico<3;ico++)gc2[ico]=gc2f[ico];
587             Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
588             if(phi2<0)phi2=2*TMath::Pi()+phi2;
589             Double_t diff = TMath::Abs(phi2-phi1); 
590             if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; 
591             if(optCuts==0 && diff<fDiffPhiforPileup){
592               Double_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
593               Double_t zc1=gc1[2];
594               Double_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
595               Double_t zc2=gc2[2];
596               Double_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius
597               fZHisto->Fill(zr0);
598             }
599             if(diff>deltaPhi)continue;
600             AliStrLine line(gc1,gc2,kTRUE);
601             Double_t cp[3];
602             Int_t retcode = line.Cross(&zeta,cp);
603             if(retcode<0)continue;
604             Double_t dca = line.GetDCA(&zeta);
605             if(dca<0.) continue;
606             if(dca>deltaR)continue;
607             Double_t deltaZ=cp[2]-zvert;
608             if(TMath::Abs(deltaZ)>dZmax)continue;
609
610             if(nolines == 0){
611               if(fLines.GetEntriesFast()>0)fLines.Clear("C");
612             }
613             Float_t cov[6];
614             recp2->GetGlobalCov(cov);
615
616
617             Double_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
618             Double_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
619             Double_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction 
620
621             Double_t curvErr=0;
622             if(bField>0.00001){
623               Double_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm 
624               Double_t dRad=TMath::Sqrt((gc1[0]-gc2[0])*(gc1[0]-gc2[0])+(gc1[1]-gc2[1])*(gc1[1]-gc2[1]));
625               Double_t aux=dRad/2.+rad1;
626               curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm
627             }
628             Double_t sigmasq[3];
629             sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor;
630             sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor;
631             sigmasq[2]=cov[5]*factor*factor;
632
633             // Multiple scattering
634             Double_t pOverMass=fMeanPSelTrk/0.140;
635             Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass);
636             Double_t p2=fMeanPSelTrk*fMeanPSelTrk;
637             Double_t rBP=GetPipeRadius();
638             Double_t dBP=0.08/35.3; // 800 um of Be
639             Double_t dL1=0.01; //approx. 1% of radiation length  
640             Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP;
641             Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1;
642             Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1));
643             Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP));
644             for(Int_t ico=0; ico<3;ico++){    
645               sigmasq[ico]+=rtantheta1*rtantheta1*factor*factor/3.;
646               sigmasq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.;
647             }
648             Double_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.};
649             if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0];
650             if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1];
651             if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2];
652             new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2);
653
654           }
655           fDetTypeRec->ResetRecPoints();
656         }
657       }
658     }
659     prpl1.Clear();
660   }
661   if(nolines == 0)return -2;
662   return nolines;
663 }
664
665 //______________________________________________________________________
666 Int_t  AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
667   // Finds the 3D vertex information using tracklets
668   Int_t retcode = -1;
669
670   Double_t xbeam=GetNominalPos()[0];
671   Double_t ybeam=GetNominalPos()[1];
672   Double_t zvert=0.;
673   Double_t deltaR=fCoarseMaxRCut;
674   Double_t dZmax=fZCutDiamond;
675   if(optCuts==1){
676     xbeam=fVert3D.GetXv();
677     ybeam=fVert3D.GetYv();
678     zvert=fVert3D.GetZv();
679     deltaR=fMaxRCut;
680     dZmax=fMaxZCut;
681     if(fPileupAlgo == 2){ 
682       dZmax=fZCutDiamond;
683       deltaR=fMaxRCut2;
684     }
685   }else if(optCuts==2){
686     xbeam=fVert3D.GetXv();
687     ybeam=fVert3D.GetYv();
688     deltaR=fMaxRCut;
689   }
690
691   Double_t rl=-fCoarseMaxRCut;
692   Double_t rh=fCoarseMaxRCut;
693   Double_t zl=-fZCutDiamond;
694   Double_t zh=fZCutDiamond;
695   Int_t nbr=(Int_t)((rh-rl)/fBinSizeR+0.0001);
696   Int_t nbz=(Int_t)((zh-zl)/fBinSizeZ+0.0001);
697   Int_t nbrcs=(Int_t)((rh-rl)/(fBinSizeR*2.)+0.0001);
698   Int_t nbzcs=(Int_t)((zh-zl)/(fBinSizeZ*2.)+0.0001);
699
700   TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
701   TH3F *h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh);
702
703   // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
704   Int_t *validate = new Int_t [fLines.GetEntriesFast()];
705   for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0;
706   for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){
707     AliStrLine *l1 = (AliStrLine*)fLines.At(i);
708     for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){
709       AliStrLine *l2 = (AliStrLine*)fLines.At(j);
710       Double_t dca=l1->GetDCA(l2);
711       if(dca > fDCAcut || dca<0.00001) continue;
712       Double_t point[3];
713       Int_t retc = l1->Cross(l2,point);
714       if(retc<0)continue;
715       Double_t deltaZ=point[2]-zvert;
716      if(TMath::Abs(deltaZ)>dZmax)continue;
717       Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
718       if(rad>fCoarseMaxRCut)continue;
719       Double_t deltaX=point[0]-xbeam;
720       Double_t deltaY=point[1]-ybeam;
721       Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
722       if(raddist>deltaR)continue;
723       validate[i]=1;
724       validate[j]=1;
725       h3d->Fill(point[0],point[1],point[2]);
726       h3dcs->Fill(point[0],point[1],point[2]);
727     }
728   }
729
730
731
732   Int_t numbtracklets=0;
733   for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++;
734   if(numbtracklets<2){
735     delete [] validate; 
736     delete h3d; 
737     delete h3dcs; 
738     return retcode; 
739   }
740
741   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
742     if(validate[i]<1)fLines.RemoveAt(i);
743   }
744   fLines.Compress();
745   AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast()));
746   delete [] validate;
747
748   // Exit here if Pileup Algorithm 2 has been chosen during second loop
749   if(fPileupAlgo == 2 && optCuts==1){
750     delete h3d; 
751     delete h3dcs;     
752     return 0;
753   }
754
755   //        Find peaks in histos
756
757   Double_t peak[3]={0.,0.,0.};
758   Int_t ntrkl,ntimes;
759   FindPeaks(h3d,peak,ntrkl,ntimes);  
760   delete h3d;
761   Double_t binsizer=(rh-rl)/nbr;
762   Double_t binsizez=(zh-zl)/nbz;
763   if(optCuts==0 && (ntrkl<=2 || ntimes>1)){
764     ntrkl=0;
765     ntimes=0;
766     FindPeaks(h3dcs,peak,ntrkl,ntimes);  
767     binsizer=(rh-rl)/nbrcs;
768     binsizez=(zh-zl)/nbzcs;
769     if(ntrkl==1 || ntimes>1){delete h3dcs; return retcode;}
770   }
771   delete h3dcs;
772
773   //         Second selection loop
774
775   Double_t bs=(binsizer+binsizez)/2.;
776   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
777     AliStrLine *l1 = (AliStrLine*)fLines.At(i);
778     if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i);
779   }
780   fLines.Compress();
781   AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast()));
782
783   if(fLines.GetEntriesFast()>1){
784     retcode=0;
785     //  find a first candidate for the primary vertex
786     fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); 
787     // make a further selection on tracklets based on this first candidate
788     fVert3D.GetXYZ(peak);
789     AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
790     for(Int_t i=0; i<fLines.GetEntriesFast();i++){
791       AliStrLine *l1 = (AliStrLine*)fLines.At(i);
792       if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i);
793     }
794     fLines.Compress();
795     AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast()));
796     if(fLines.GetEntriesFast()>1){// this new tracklet selection is used
797       fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
798     }
799   }
800   return retcode;  
801 }
802
803 //________________________________________________________
804 void AliITSVertexer3D::SetMeanPPtSelTracks(Double_t fieldTesla){
805   // Sets mean values of Pt based on the field
806   // for P (used in multiple scattering) the most probable value is used
807   if(TMath::Abs(fieldTesla-0.5)<0.01){
808     SetMeanPSelTracks(0.375);
809     SetMeanPtSelTracks(0.630);
810   }else if(TMath::Abs(fieldTesla-0.4)<0.01){
811     SetMeanPSelTracks(0.375);
812     SetMeanPtSelTracks(0.580);
813   }else if(TMath::Abs(fieldTesla-0.2)<0.01){
814     SetMeanPSelTracks(0.375);
815     SetMeanPtSelTracks(0.530);
816   }else if(fieldTesla<0.00001){
817     SetMeanPSelTracks(0.375);
818     SetMeanPtSelTracks(0.230);
819   }else{
820     SetMeanPSelTracks();
821     SetMeanPtSelTracks();
822   }
823 }
824
825 //________________________________________________________
826 void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){
827   // Finds bin with max contents in 3D histo of tracket intersections
828   TAxis *xax = histo->GetXaxis();  
829   TAxis *yax = histo->GetYaxis();
830   TAxis *zax = histo->GetZaxis();
831   peak[0]=0.;
832   peak[1]=0.;
833   peak[2]=0.;
834   nOfTracklets = 0;
835   nOfTimes=0;
836   Int_t peakbin[3]={0,0,0};
837   Int_t peak2bin[3]={-1,-1,-1};
838   Int_t bc2=-1;
839   for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
840     Double_t xval = xax->GetBinCenter(i);
841     for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
842       Double_t yval = yax->GetBinCenter(j);
843       for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
844         Double_t zval = zax->GetBinCenter(k);
845         Int_t bc =(Int_t)histo->GetBinContent(i,j,k);
846         if(bc==0) continue;
847         if(bc>nOfTracklets){
848           nOfTracklets=bc;
849           peak[2] = zval;
850           peak[1] = yval;
851           peak[0] = xval;
852           peakbin[2] = k;
853           peakbin[1] = j;
854           peakbin[0] = i;
855           peak2bin[2] = -1;
856           peak2bin[1] = -1;
857           peak2bin[0] = -1;
858           bc2=-1;
859           nOfTimes = 1;
860         }else if(bc==nOfTracklets){
861           if(TMath::Abs(i-peakbin[0])<=1 && TMath::Abs(j-peakbin[1])<=1 && TMath::Abs(k-peakbin[2])<=1){
862             peak2bin[2] = k;
863             peak2bin[1] = j;
864             peak2bin[0] = i;
865             bc2=bc;
866             nOfTimes = 1;
867           }else{
868             nOfTimes++;
869           }
870         }
871       }
872     }
873   }
874   if(peak2bin[0]>=-1 && bc2!=-1){ // two contiguous peak-cells with same contents
875     peak[0]=0.5*(xax->GetBinCenter(peakbin[0])+xax->GetBinCenter(peak2bin[0]));
876     peak[1]=0.5*(yax->GetBinCenter(peakbin[1])+yax->GetBinCenter(peak2bin[1]));
877     peak[2]=0.5*(zax->GetBinCenter(peakbin[2])+zax->GetBinCenter(peak2bin[2]));
878     nOfTracklets+=bc2;
879     nOfTimes=1;
880   }
881 }
882 //________________________________________________________
883 void AliITSVertexer3D::MarkUsedClusters(){
884   // Mark clusters of tracklets used in vertex claulation
885   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
886     AliStrLine *lin = (AliStrLine*)fLines.At(i);
887     Int_t idClu1=lin->GetIdPoint(0);
888     Int_t idClu2=lin->GetIdPoint(1);
889     fUsedCluster.SetBitNumber(idClu1);
890     fUsedCluster.SetBitNumber(idClu2);
891   }
892 }
893 //________________________________________________________
894 Int_t AliITSVertexer3D::RemoveTracklets(){
895   // Remove trackelts close to first found vertex
896   Double_t vert[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
897   Int_t nRemoved=0;
898   for(Int_t i=0; i<fLines.GetEntriesFast();i++){
899     AliStrLine *lin = (AliStrLine*)fLines.At(i);
900     if(lin->GetDistFromPoint(vert)<fDCAforPileup){
901       Int_t idClu1=lin->GetIdPoint(0);
902       Int_t idClu2=lin->GetIdPoint(1);
903       fUsedCluster.SetBitNumber(idClu1);
904       fUsedCluster.SetBitNumber(idClu2);
905       fLines.RemoveAt(i);
906       ++nRemoved;
907     }
908   }
909   fLines.Compress();
910   return nRemoved;
911 }
912 //________________________________________________________
913 void AliITSVertexer3D::FindOther3DVertices(TTree *itsClusterTree){
914   // pileup identification based on 3D vertexing with not used clusters
915   MarkUsedClusters();
916   fLines.Clear("C");
917   Int_t nolines = FindTracklets(itsClusterTree,2);
918   if(nolines>=2){
919     Int_t nr=RemoveTracklets();
920     nolines-=nr;
921     if(nolines>=2){
922       Int_t rc=Prepare3DVertex(2);
923       if(rc==0){ 
924         fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
925         if(fVert3D.GetNContributors()>=fMinTrackletsForPilup){
926           fIsPileup=kTRUE;
927           fNoVertices=2;
928           fVertArray = new AliESDVertex[2];
929           fVertArray[0]=(*fCurrentVertex);
930           fVertArray[1]=fVert3D;
931           fZpuv=fVert3D.GetZv();
932           fNTrpuv=fVert3D.GetNContributors();
933         }
934       }
935     }
936   }
937 }
938 //______________________________________________________________________
939 void AliITSVertexer3D::PileupFromZ(){
940   // Calls the pileup algorithm of ALiITSVertexerZ
941   Int_t binmin, binmax;
942   Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax);   
943   if(nPeaks==2)AliWarning("2 peaks found");
944   Int_t firstPeakCont=0;
945   Double_t firstPeakPos=0.;
946   for(Int_t i=binmin-1;i<=binmax+1;i++){
947     firstPeakCont+=static_cast<Int_t>(fZHisto->GetBinContent(i));
948     firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i);
949   }
950   if(firstPeakCont>0){ 
951     firstPeakPos/=firstPeakCont;
952     Int_t ncontr2=0;
953     if(firstPeakCont>fMinTrackletsForPilup){     
954       Float_t secPeakPos;
955       ncontr2=AliITSVertexerZ::FindSecondPeak(fZHisto,binmin,binmax,secPeakPos);
956       if(ncontr2>=fMinTrackletsForPilup){ 
957         fIsPileup=kTRUE;
958         fNoVertices=2;
959         AliESDVertex secondVert(secPeakPos,0.1,ncontr2);
960         fVertArray = new AliESDVertex[2];
961         fVertArray[0]=(*fCurrentVertex);
962         fVertArray[1]=secondVert;
963         fZpuv=secPeakPos;
964         fNTrpuv=ncontr2;
965       }
966     }
967   }
968 }
969 //________________________________________________________
970 void AliITSVertexer3D::PrintStatus() const {
971   // Print current status
972   printf("========= First step selections =====================\n");
973   printf("Cut on diamond (Z) %f\n",fZCutDiamond);
974   printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut);
975   printf("Loose cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut);
976   printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut);
977   printf("========= Second step selections ====================\n");
978   printf("Cut on tracklet-to-first-vertex Z distance %f\n",fMaxZCut);
979   printf("Max Phi difference: %f\n",fDiffPhiMax);
980   printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut);
981   printf("Cut on tracklet DCA to beam axis (algo2) %f\n",fMaxRCut2);
982   printf("========= Pileup selections =========================\n");
983   printf("Pileup algo: %d\n",fPileupAlgo);
984   printf("Min DCA to 1st vertex for pileup (algo 0 and 1): %f\n",fDCAforPileup);
985   printf("Cut on distance between pair-vertices  (algo 2): %f\n",fCutOnPairs);
986   printf("=======================================================\n");
987 }