[u/mrichter/AliRoot.git] / ITS / AliITSVertexerCosmics.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#include <TClonesArray.h>
#include "AliLog.h"
#include "AliESDVertex.h"
#include "AliRunLoader.h"
#include "AliITSLoader.h"
#include "AliITSgeomTGeo.h"
#include "AliITSRecPoint.h"
#include "AliITSReconstructor.h"
#include "AliITSVertexerCosmics.h"
#include "AliStrLine.h"

//------------------------------------------------------------------------
// This class implements a method to construct a "fake" primary
// vertex for cosmic events in which the muon crosses one of 5 inner
// ITS layers. A fake primary vertex is needed for the reconstruction,
// with e.g. AliITStrackerSA, of the two tracks produced by the muon 
// in the ITS.
//   We build pairs of clusters on a given layer and define the fake vertex as
// the mid-point of the straight line joining the two clusters.
//   We use the innermost layer that has at least two clusters.
//   We reject the background by requiring at least one cluster on the outer
// layer, closer than fMaxDistOnOuterLayer to the tracklet prolongation.
//   We can reject (potentially pathological) events with the muon track
// tangential to the layer by the requiring the radial position of
// the vertex to be smaller than fMaxVtxRadius.
//   Due to background clusters, more than one vertex per event can 
// be found. We consider the first found.
//   The errors on x,y,z of the vertex are calculated as errors on the mean
// of clusters coordinates. Non-diag elements of vertex cov. mat. are set to 0.
//   The number of contributors set in the AliESDVertex object is the
// number of the layer on which the tracklet was built; if this number is -1, 
// the procedure could not find a vertex position and by default 
// the vertex coordinates are set to (0,0,0) with large errors (100,100,100)
//
// Origin: A.Dainese, andrea.dainese@lnl.infn.it
//-------------------------------------------------------------------------

ClassImp(AliITSVertexerCosmics)

//-------------------------------------------------------------------------
AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
fMaxDistOnOuterLayer(0),
fMinDist2Vtxs(0)
{
  // Default constructor
  SetFirstLastModules(0,0,79);
  SetFirstLastModules(1,80,239);
  SetFirstLastModules(2,240,323);
  SetFirstLastModules(3,324,499);
  SetFirstLastModules(4,500,1247);
  SetFirstLastModules(5,1248,2197);
  /*
  SetMaxVtxRadius(0,3.5);
  SetMaxVtxRadius(1,6.5);
  SetMaxVtxRadius(2,14.5);
  SetMaxVtxRadius(3,23.5);
  SetMaxVtxRadius(4,37.5);
  SetMaxVtxRadius(5,42.5);
  */  
  SetMaxVtxRadius(0,5.5);
  SetMaxVtxRadius(1,8.5);
  SetMaxVtxRadius(2,18.5);
  SetMaxVtxRadius(3,28.5);
  SetMaxVtxRadius(4,39.5);
  SetMaxVtxRadius(5,48.5);
  
  SetMaxDistOnOuterLayer();
  SetMinDist2Vtxs();
}
//--------------------------------------------------------------------------
AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(Int_t evnumber) 
{
  // Defines the AliESDVertex for the current event

  fCurrentVertex = 0;
  AliRunLoader *rl =AliRunLoader::GetRunLoader();
  AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
  itsLoader->LoadRecPoints();
  rl->GetEvent(evnumber);

  TTree *rpTree = itsLoader->TreeR();

  TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
  rpTree->SetBranchAddress("ITSRecPoints",&recpoints);

  Int_t lay,lad,det; 

  Double_t pos[3]={0.,0.,0.};
  Double_t err[3]={100.,100.,100.};
  Int_t ncontributors = -1;

  // Search for innermost layer with at least two clusters 
  // on two different modules
  Int_t ilayer=0,ilayer2=0;
  while(ilayer<6) {
    if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) {
      ilayer++;
      continue;
    }
    Int_t nHitModules=0;
    for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer]; imodule++) {
      rpTree->GetEvent(imodule);
      AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
      lay -= 1;  // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
      if(lay!=ilayer) AliFatal("Layer mismatch!");
      if(recpoints->GetEntriesFast()>0) nHitModules++;
    }
    if(nHitModules>=2) break;
    ilayer++;
  }

  ilayer2=ilayer+1;
  while(ilayer2<6) {
    if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer2)) break;
    ilayer2++;
  }

  // try tracklet on SPD2 and point on SPD1
  if(ilayer==1 && ilayer2>5 && !AliITSReconstructor::GetRecoParam()->GetLayersToSkip(0)) {ilayer=0; ilayer2=1;}

  if(ilayer>4 || ilayer2>5) {
    AliWarning("Not enough clusters");
    delete recpoints;
    itsLoader->UnloadRecPoints();
    fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
    fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
    fCurrentVertex->SetNContributors(ncontributors);
    return fCurrentVertex;
  }


  const Int_t arrSize = 1000;
  Float_t xclInnLay[arrSize],yclInnLay[arrSize],zclInnLay[arrSize],modclInnLay[arrSize];
  Float_t e2xclInnLay[arrSize],e2yclInnLay[arrSize],e2zclInnLay[arrSize];
  Float_t e2xclOutLay[arrSize],e2yclOutLay[arrSize],e2zclOutLay[arrSize];
  Int_t nclInnLayStored=0;
  Float_t xclOutLay[arrSize],yclOutLay[arrSize],zclOutLay[arrSize],modclOutLay[arrSize];
  Int_t nclOutLayStored=0;
  Int_t nRecPoints,nRecPointsInnLay=0;

  Float_t gc[3],gcov[5];

  Float_t matchOutLayValue;
  Float_t distxyInnLay,distxyInnLayBest=0.;
  Double_t p1[3],p2[3],p3[3];

  Float_t xvtx,yvtx,zvtx,rvtx;
  Int_t i1InnLayBest=-1,i2InnLayBest=-1;


  // Collect clusters in the selected layer and the outer one
  for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer2]; imodule++) {
    rpTree->GetEvent(imodule);
    AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
    lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
    nRecPoints=recpoints->GetEntriesFast();
    if(imodule<=fLast[ilayer]) nRecPointsInnLay += nRecPoints;
    //printf("cosmics: module %d clusters %d\n",imodule,nRecPoints);
    for(Int_t irp=0; irp<nRecPoints; irp++) {
      AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
      // Local coordinates of this recpoint
      rp->GetGlobalXYZ(gc);
      if(lay==ilayer) { // store InnLay clusters
	xclInnLay[nclInnLayStored]=gc[0];
	yclInnLay[nclInnLayStored]=gc[1];
	zclInnLay[nclInnLayStored]=gc[2];
	rp->GetGlobalCov(gcov);
	e2xclInnLay[nclInnLayStored]=gcov[0];
	e2yclInnLay[nclInnLayStored]=gcov[3];
	e2zclInnLay[nclInnLayStored]=gcov[5];
	modclInnLay[nclInnLayStored]=imodule;
	nclInnLayStored++;
      }
      if(lay==ilayer2) { // store OutLay clusters
	xclOutLay[nclOutLayStored]=gc[0];
	yclOutLay[nclOutLayStored]=gc[1];
	zclOutLay[nclOutLayStored]=gc[2];
	rp->GetGlobalCov(gcov);
	e2xclOutLay[nclOutLayStored]=gcov[0];
	e2yclOutLay[nclOutLayStored]=gcov[3];
	e2zclOutLay[nclOutLayStored]=gcov[5];
	modclOutLay[nclOutLayStored]=imodule;
	nclOutLayStored++;
      }
      if(nclInnLayStored>arrSize || nclOutLayStored>arrSize) {
	//AliFatal("More than arrSize clusters per layer");
	AliWarning("Too many clusters per layer");
	delete recpoints;
	itsLoader->UnloadRecPoints();
	fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
	fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
	fCurrentVertex->SetNContributors(ncontributors);
	return fCurrentVertex;
      }
    }// end clusters in a module
  }// end modules


  Int_t i1InnLay,i2InnLay,iOutLay;

  // build fake vertices
  printf("Building tracklets on layer %d\n",ilayer);

  // InnLay - first cluster
  for(i1InnLay=0; i1InnLay<nclInnLayStored; i1InnLay++) { 
    p1[0]=xclInnLay[i1InnLay]; 
    p1[1]=yclInnLay[i1InnLay]; 
    p1[2]=zclInnLay[i1InnLay];
    // InnLay - second cluster
    for(i2InnLay=i1InnLay+1; i2InnLay<nclInnLayStored; i2InnLay++) { 
      if(modclInnLay[i1InnLay]==modclInnLay[i2InnLay]) continue;
      p2[0]=xclInnLay[i2InnLay]; 
      p2[1]=yclInnLay[i2InnLay]; 
      p2[2]=zclInnLay[i2InnLay];
      // look for point on OutLay
      AliStrLine innLayline(p1,p2,kTRUE);
      for(iOutLay=0; iOutLay<nclOutLayStored; iOutLay++) {
	p3[0]=xclOutLay[iOutLay]; 
	p3[1]=yclOutLay[iOutLay]; 
	p3[2]=zclOutLay[iOutLay];
	//printf("(%f,%f) (%f,%f)     (%f,%f) %f\n",p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],innLayline.GetDistFromPoint(p3));
	matchOutLayValue=innLayline.GetDistFromPoint(p3);
	distxyInnLay = (p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1]);
	if(matchOutLayValue<fMaxDistOnOuterLayer &&
	   distxyInnLay>distxyInnLayBest) { 
	  //printf("found\n");
	  distxyInnLayBest=distxyInnLay;
	  i1InnLayBest=i1InnLay;
	  i2InnLayBest=i2InnLay;
	}
      }
    } // InnLay - second cluster
  } // InnLay - first cluster

  if(i1InnLayBest>-1 && i2InnLayBest>-1) { 
    xvtx = 0.5*(xclInnLay[i1InnLayBest]+xclInnLay[i2InnLayBest]);
    yvtx = 0.5*(yclInnLay[i1InnLayBest]+yclInnLay[i2InnLayBest]);
    zvtx = 0.5*(zclInnLay[i1InnLayBest]+zclInnLay[i2InnLayBest]);
    rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
    if(rvtx<fMaxVtxRadius[ilayer]) {
      ncontributors = ilayer;
      pos[0] = xvtx;
      pos[1] = yvtx;
      pos[2] = zvtx;
      err[0]=TMath::Sqrt(0.25*(e2xclInnLay[i1InnLayBest]+e2xclInnLay[i2InnLayBest])); 
      err[1]=TMath::Sqrt(0.25*(e2yclInnLay[i1InnLayBest]+e2yclInnLay[i2InnLayBest])); 
      err[2]=TMath::Sqrt(0.25*(e2zclInnLay[i1InnLayBest]+e2zclInnLay[i2InnLayBest]));
    }

  } else { // give it a try exchanging InnLay and OutLay

    // OutLay - first cluster
    for(i1InnLay=0; i1InnLay<nclOutLayStored; i1InnLay++) { 
      p1[0]=xclOutLay[i1InnLay]; 
      p1[1]=yclOutLay[i1InnLay]; 
      p1[2]=zclOutLay[i1InnLay];
      // OutLay - second cluster
      for(i2InnLay=i1InnLay+1; i2InnLay<nclOutLayStored; i2InnLay++) { 
	if(modclOutLay[i1InnLay]==modclOutLay[i2InnLay]) continue;
	p2[0]=xclOutLay[i2InnLay]; 
	p2[1]=yclOutLay[i2InnLay]; 
	p2[2]=zclOutLay[i2InnLay];
	// look for point on InnLay
	AliStrLine outLayline(p1,p2,kTRUE);
	for(iOutLay=0; iOutLay<nclInnLayStored; iOutLay++) {
	  p3[0]=xclInnLay[iOutLay]; 
	  p3[1]=yclInnLay[iOutLay]; 
	  p3[2]=zclInnLay[iOutLay];
	  //printf(" %f\n",InnLayline.GetDistFromPoint(p3));
	  matchOutLayValue=outLayline.GetDistFromPoint(p3);
	  distxyInnLay = (p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1]);
	  if(matchOutLayValue<fMaxDistOnOuterLayer &&
	     distxyInnLay>distxyInnLayBest) { 
	    distxyInnLayBest=distxyInnLay;
	    i1InnLayBest=i1InnLay;
	    i2InnLayBest=i2InnLay;
	  }
	}
      } // OutLay - second cluster
    } // OutLay - first cluster

    if(i1InnLayBest>-1 && i2InnLayBest>-1) { 
      xvtx = 0.5*(xclOutLay[i1InnLayBest]+xclOutLay[i2InnLayBest]);
      yvtx = 0.5*(yclOutLay[i1InnLayBest]+yclOutLay[i2InnLayBest]);
      zvtx = 0.5*(zclOutLay[i1InnLayBest]+zclOutLay[i2InnLayBest]);
      rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
      if(rvtx<fMaxVtxRadius[ilayer]) {
	ncontributors = ilayer2;
	pos[0] = xvtx;
	pos[1] = yvtx;
	pos[2] = zvtx;
	err[0]=TMath::Sqrt(0.25*(e2xclOutLay[i1InnLayBest]+e2xclOutLay[i2InnLayBest])); 
	err[1]=TMath::Sqrt(0.25*(e2yclOutLay[i1InnLayBest]+e2yclOutLay[i2InnLayBest])); 
	err[2]=TMath::Sqrt(0.25*(e2zclOutLay[i1InnLayBest]+e2zclOutLay[i2InnLayBest]));
      }
    }
  } // give it a try exchanging InnLay and OutLay

  fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
  fCurrentVertex->SetTitle("cosmics fake vertex");
  fCurrentVertex->SetNContributors(ncontributors);
  //fCurrentVertex->Print();

  delete recpoints;
  itsLoader->UnloadRecPoints();

  return fCurrentVertex;
}  
//-------------------------------------------------------------------------
void AliITSVertexerCosmics::FindVertices()
{
  // computes the vertices of the events in the range FirstEvent - LastEvent
  AliRunLoader *rl = AliRunLoader::GetRunLoader();
  AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  itsLoader->ReloadRecPoints();
  for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
    //  printf("Processing event %d\n",i);
    rl->GetEvent(i);
    FindVertexForCurrentEvent(i);
    if(fCurrentVertex){
      WriteCurrentVertex();
    }
  }
}
//-------------------------------------------------------------------------
void AliITSVertexerCosmics::PrintStatus() const 
{
  // Print current status
  printf("=======================================================\n");
  printf(" fMaxDistOnOuterLayer: %f\n",fMaxDistOnOuterLayer);
  printf(" fMaxVtxRadius[0]:  %f\n",fMaxVtxRadius[0]);
  printf(" fMinDist2Vtxs:  %f\n",fMinDist2Vtxs);
  printf("=======================================================\n");
}
//-------------------------------------------------------------------------
Commit	Line	Data
3acc14d5	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, 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	#include <TClonesArray.h>
	17	#include "AliLog.h"
	18	#include "AliESDVertex.h"
	19	#include "AliRunLoader.h"
	20	#include "AliITSLoader.h"
0333e459	21	#include "AliITSgeomTGeo.h"
3acc14d5	22	#include "AliITSRecPoint.h"
81ba01de	23	#include "AliITSReconstructor.h"
3acc14d5	24	#include "AliITSVertexerCosmics.h"
	25	#include "AliStrLine.h"
	26
	27	//------------------------------------------------------------------------
	28	// This class implements a method to construct a "fake" primary
b8ed1a92	29	// vertex for cosmic events in which the muon crosses one of 5 inner
b8ed1a92	30	// ITS layers. A fake primary vertex is needed for the reconstruction,
3acc14d5	31	// with e.g. AliITStrackerSA, of the two tracks produced by the muon
3acc14d5	32	// in the ITS.
b8ed1a92	33	// We build pairs of clusters on a given layer and define the fake vertex as
3acc14d5	34	// the mid-point of the straight line joining the two clusters.
b8ed1a92	35	// We use the innermost layer that has at least two clusters.
	36	// We reject the background by requiring at least one cluster on the outer
	37	// layer, closer than fMaxDistOnOuterLayer to the tracklet prolongation.
3acc14d5	38	// We can reject (potentially pathological) events with the muon track
b8ed1a92	39	// tangential to the layer by the requiring the radial position of
3acc14d5	40	// the vertex to be smaller than fMaxVtxRadius.
	41	// Due to background clusters, more than one vertex per event can
	42	// be found. We consider the first found.
b8ed1a92	43	// The errors on x,y,z of the vertex are calculated as errors on the mean
b8ed1a92	44	// of clusters coordinates. Non-diag elements of vertex cov. mat. are set to 0.
3acc14d5	45	// The number of contributors set in the AliESDVertex object is the
b8ed1a92	46	// number of the layer on which the tracklet was built; if this number is -1,
3acc14d5	47	// the procedure could not find a vertex position and by default
3acc14d5	48	// the vertex coordinates are set to (0,0,0) with large errors (100,100,100)
3acc14d5	49	//
	50	// Origin: A.Dainese, andrea.dainese@lnl.infn.it
	51	//-------------------------------------------------------------------------
	52
	53	ClassImp(AliITSVertexerCosmics)
	54
	55	//-------------------------------------------------------------------------
	56	AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
b8ed1a92	57	fMaxDistOnOuterLayer(0),
3acc14d5	58	fMinDist2Vtxs(0)
	59	{
	60	// Default constructor
b8ed1a92	61	SetFirstLastModules(0,0,79);
	62	SetFirstLastModules(1,80,239);
	63	SetFirstLastModules(2,240,323);
	64	SetFirstLastModules(3,324,499);
	65	SetFirstLastModules(4,500,1247);
	66	SetFirstLastModules(5,1248,2197);
7e33c4e6	67	/*
b8ed1a92	68	SetMaxVtxRadius(0,3.5);
	69	SetMaxVtxRadius(1,6.5);
	70	SetMaxVtxRadius(2,14.5);
	71	SetMaxVtxRadius(3,23.5);
	72	SetMaxVtxRadius(4,37.5);
	73	SetMaxVtxRadius(5,42.5);
7e33c4e6	74	*/
81ba01de	75	SetMaxVtxRadius(0,5.5);
	76	SetMaxVtxRadius(1,8.5);
	77	SetMaxVtxRadius(2,18.5);
	78	SetMaxVtxRadius(3,28.5);
	79	SetMaxVtxRadius(4,39.5);
	80	SetMaxVtxRadius(5,48.5);
7e33c4e6	81
b8ed1a92	82	SetMaxDistOnOuterLayer();
3acc14d5	83	SetMinDist2Vtxs();
	84	}
	85	//--------------------------------------------------------------------------
	86	AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(Int_t evnumber)
	87	{
	88	// Defines the AliESDVertex for the current event
	89
	90	fCurrentVertex = 0;
	91	AliRunLoader *rl =AliRunLoader::GetRunLoader();
	92	AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
3acc14d5	93	itsLoader->LoadRecPoints();
	94	rl->GetEvent(evnumber);
	95
	96	TTree *rpTree = itsLoader->TreeR();
	97
	98	TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
	99	rpTree->SetBranchAddress("ITSRecPoints",&recpoints);
	100
3acc14d5	101	Int_t lay,lad,det;
3acc14d5	102
81ba01de	103	Double_t pos[3]={0.,0.,0.};
	104	Double_t err[3]={100.,100.,100.};
	105	Int_t ncontributors = -1;
	106
b8ed1a92	107	// Search for innermost layer with at least two clusters
b8ed1a92	108	// on two different modules
81ba01de	109	Int_t ilayer=0,ilayer2=0;
b8ed1a92	110	while(ilayer<6) {
81ba01de	111	if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) {
	112	ilayer++;
	113	continue;
	114	}
b8ed1a92	115	Int_t nHitModules=0;
	116	for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer]; imodule++) {
	117	rpTree->GetEvent(imodule);
	118	AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
	119	lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
	120	if(lay!=ilayer) AliFatal("Layer mismatch!");
	121	if(recpoints->GetEntriesFast()>0) nHitModules++;
	122	}
	123	if(nHitModules>=2) break;
	124	ilayer++;
	125	}
81ba01de	126
	127	ilayer2=ilayer+1;
	128	while(ilayer2<6) {
	129	if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer2)) break;
	130	ilayer2++;
	131	}
	132
7e33c4e6	133	// try tracklet on SPD2 and point on SPD1
7e33c4e6	134	if(ilayer==1 && ilayer2>5 && !AliITSReconstructor::GetRecoParam()->GetLayersToSkip(0)) {ilayer=0; ilayer2=1;}
81ba01de	135
	136	if(ilayer>4 \|\| ilayer2>5) {
	137	AliWarning("Not enough clusters");
	138	delete recpoints;
	139	itsLoader->UnloadRecPoints();
	140	fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
	141	fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
	142	fCurrentVertex->SetNContributors(ncontributors);
	143	return fCurrentVertex;
	144	}
	145
b8ed1a92	146
a1fca0e6	147	const Int_t arrSize = 1000;
	148	Float_t xclInnLay[arrSize],yclInnLay[arrSize],zclInnLay[arrSize],modclInnLay[arrSize];
	149	Float_t e2xclInnLay[arrSize],e2yclInnLay[arrSize],e2zclInnLay[arrSize];
81ba01de	150	Float_t e2xclOutLay[arrSize],e2yclOutLay[arrSize],e2zclOutLay[arrSize];
b8ed1a92	151	Int_t nclInnLayStored=0;
a1fca0e6	152	Float_t xclOutLay[arrSize],yclOutLay[arrSize],zclOutLay[arrSize],modclOutLay[arrSize];
b8ed1a92	153	Int_t nclOutLayStored=0;
	154	Int_t nRecPoints,nRecPointsInnLay=0;
	155
	156	Float_t gc[3],gcov[5];
	157
81ba01de	158	Float_t matchOutLayValue;
81ba01de	159	Float_t distxyInnLay,distxyInnLayBest=0.;
b8ed1a92	160	Double_t p1[3],p2[3],p3[3];
81ba01de	161
b8ed1a92	162	Float_t xvtx,yvtx,zvtx,rvtx;
81ba01de	163	Int_t i1InnLayBest=-1,i2InnLayBest=-1;
3acc14d5	164
b107ee21	165
b8ed1a92	166	// Collect clusters in the selected layer and the outer one
81ba01de	167	for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer2]; imodule++) {
3acc14d5	168	rpTree->GetEvent(imodule);
0333e459	169	AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
b8ed1a92	170	lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
	171	nRecPoints=recpoints->GetEntriesFast();
	172	if(imodule<=fLast[ilayer]) nRecPointsInnLay += nRecPoints;
	173	//printf("cosmics: module %d clusters %d\n",imodule,nRecPoints);
	174	for(Int_t irp=0; irp<nRecPoints; irp++) {
3acc14d5	175	AliITSRecPoint rp=(AliITSRecPoint)recpoints->UncheckedAt(irp);
3acc14d5	176	// Local coordinates of this recpoint
b8ed1a92	177	rp->GetGlobalXYZ(gc);
	178	if(lay==ilayer) { // store InnLay clusters
	179	xclInnLay[nclInnLayStored]=gc[0];
	180	yclInnLay[nclInnLayStored]=gc[1];
	181	zclInnLay[nclInnLayStored]=gc[2];
	182	rp->GetGlobalCov(gcov);
	183	e2xclInnLay[nclInnLayStored]=gcov[0];
	184	e2yclInnLay[nclInnLayStored]=gcov[3];
	185	e2zclInnLay[nclInnLayStored]=gcov[5];
	186	modclInnLay[nclInnLayStored]=imodule;
	187	nclInnLayStored++;
3acc14d5	188	}
81ba01de	189	if(lay==ilayer2) { // store OutLay clusters
b8ed1a92	190	xclOutLay[nclOutLayStored]=gc[0];
	191	yclOutLay[nclOutLayStored]=gc[1];
	192	zclOutLay[nclOutLayStored]=gc[2];
81ba01de	193	rp->GetGlobalCov(gcov);
	194	e2xclOutLay[nclOutLayStored]=gcov[0];
	195	e2yclOutLay[nclOutLayStored]=gcov[3];
	196	e2zclOutLay[nclOutLayStored]=gcov[5];
b8ed1a92	197	modclOutLay[nclOutLayStored]=imodule;
b8ed1a92	198	nclOutLayStored++;
3acc14d5	199	}
b107ee21	200	if(nclInnLayStored>arrSize \|\| nclOutLayStored>arrSize) {
	201	//AliFatal("More than arrSize clusters per layer");
	202	AliWarning("Too many clusters per layer");
	203	delete recpoints;
	204	itsLoader->UnloadRecPoints();
	205	fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
	206	fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
	207	fCurrentVertex->SetNContributors(ncontributors);
	208	return fCurrentVertex;
	209	}
3acc14d5	210	}// end clusters in a module
b8ed1a92	211	}// end modules
3acc14d5	212
81ba01de	213
	214	Int_t i1InnLay,i2InnLay,iOutLay;
	215
3acc14d5	216	// build fake vertices
81ba01de	217	printf("Building tracklets on layer %d\n",ilayer);
81ba01de	218
b8ed1a92	219	// InnLay - first cluster
81ba01de	220	for(i1InnLay=0; i1InnLay<nclInnLayStored; i1InnLay++) {
b8ed1a92	221	p1[0]=xclInnLay[i1InnLay];
	222	p1[1]=yclInnLay[i1InnLay];
	223	p1[2]=zclInnLay[i1InnLay];
	224	// InnLay - second cluster
81ba01de	225	for(i2InnLay=i1InnLay+1; i2InnLay<nclInnLayStored; i2InnLay++) {
b8ed1a92	226	if(modclInnLay[i1InnLay]==modclInnLay[i2InnLay]) continue;
	227	p2[0]=xclInnLay[i2InnLay];
	228	p2[1]=yclInnLay[i2InnLay];
	229	p2[2]=zclInnLay[i2InnLay];
	230	// look for point on OutLay
81ba01de	231	AliStrLine innLayline(p1,p2,kTRUE);
81ba01de	232	for(iOutLay=0; iOutLay<nclOutLayStored; iOutLay++) {
b8ed1a92	233	p3[0]=xclOutLay[iOutLay];
	234	p3[1]=yclOutLay[iOutLay];
	235	p3[2]=zclOutLay[iOutLay];
7e33c4e6	236	//printf("(%f,%f) (%f,%f) (%f,%f) %f\n",p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],innLayline.GetDistFromPoint(p3));
81ba01de	237	matchOutLayValue=innLayline.GetDistFromPoint(p3);
	238	distxyInnLay = (p1[0]-p2[0])(p1[0]-p2[0])+(p1[1]-p2[1])(p1[1]-p2[1]);
	239	if(matchOutLayValue<fMaxDistOnOuterLayer &&
	240	distxyInnLay>distxyInnLayBest) {
7e33c4e6	241	//printf("found\n");
81ba01de	242	distxyInnLayBest=distxyInnLay;
	243	i1InnLayBest=i1InnLay;
	244	i2InnLayBest=i2InnLay;
	245	}
3acc14d5	246	}
b8ed1a92	247	} // InnLay - second cluster
b8ed1a92	248	} // InnLay - first cluster
3acc14d5	249
81ba01de	250	if(i1InnLayBest>-1 && i2InnLayBest>-1) {
	251	xvtx = 0.5*(xclInnLay[i1InnLayBest]+xclInnLay[i2InnLayBest]);
	252	yvtx = 0.5*(yclInnLay[i1InnLayBest]+yclInnLay[i2InnLayBest]);
	253	zvtx = 0.5*(zclInnLay[i1InnLayBest]+zclInnLay[i2InnLayBest]);
	254	rvtx = TMath::Sqrt(xvtxxvtx+yvtxyvtx);
	255	if(rvtx<fMaxVtxRadius[ilayer]) {
	256	ncontributors = ilayer;
	257	pos[0] = xvtx;
	258	pos[1] = yvtx;
	259	pos[2] = zvtx;
	260	err[0]=TMath::Sqrt(0.25*(e2xclInnLay[i1InnLayBest]+e2xclInnLay[i2InnLayBest]));
	261	err[1]=TMath::Sqrt(0.25*(e2yclInnLay[i1InnLayBest]+e2yclInnLay[i2InnLayBest]));
	262	err[2]=TMath::Sqrt(0.25*(e2zclInnLay[i1InnLayBest]+e2zclInnLay[i2InnLayBest]));
	263	}
81ba01de	264
7e33c4e6	265	} else { // give it a try exchanging InnLay and OutLay
7e33c4e6	266
81ba01de	267	// OutLay - first cluster
	268	for(i1InnLay=0; i1InnLay<nclOutLayStored; i1InnLay++) {
	269	p1[0]=xclOutLay[i1InnLay];
	270	p1[1]=yclOutLay[i1InnLay];
	271	p1[2]=zclOutLay[i1InnLay];
	272	// OutLay - second cluster
	273	for(i2InnLay=i1InnLay+1; i2InnLay<nclOutLayStored; i2InnLay++) {
	274	if(modclOutLay[i1InnLay]==modclOutLay[i2InnLay]) continue;
	275	p2[0]=xclOutLay[i2InnLay];
	276	p2[1]=yclOutLay[i2InnLay];
	277	p2[2]=zclOutLay[i2InnLay];
	278	// look for point on InnLay
	279	AliStrLine outLayline(p1,p2,kTRUE);
	280	for(iOutLay=0; iOutLay<nclInnLayStored; iOutLay++) {
	281	p3[0]=xclInnLay[iOutLay];
	282	p3[1]=yclInnLay[iOutLay];
	283	p3[2]=zclInnLay[iOutLay];
	284	//printf(" %f\n",InnLayline.GetDistFromPoint(p3));
	285	matchOutLayValue=outLayline.GetDistFromPoint(p3);
	286	distxyInnLay = (p1[0]-p2[0])(p1[0]-p2[0])+(p1[1]-p2[1])(p1[1]-p2[1]);
	287	if(matchOutLayValue<fMaxDistOnOuterLayer &&
	288	distxyInnLay>distxyInnLayBest) {
	289	distxyInnLayBest=distxyInnLay;
	290	i1InnLayBest=i1InnLay;
	291	i2InnLayBest=i2InnLay;
	292	}
	293	}
	294	} // OutLay - second cluster
	295	} // OutLay - first cluster
81ba01de	296
7e33c4e6	297	if(i1InnLayBest>-1 && i2InnLayBest>-1) {
	298	xvtx = 0.5*(xclOutLay[i1InnLayBest]+xclOutLay[i2InnLayBest]);
	299	yvtx = 0.5*(yclOutLay[i1InnLayBest]+yclOutLay[i2InnLayBest]);
	300	zvtx = 0.5*(zclOutLay[i1InnLayBest]+zclOutLay[i2InnLayBest]);
	301	rvtx = TMath::Sqrt(xvtxxvtx+yvtxyvtx);
	302	if(rvtx<fMaxVtxRadius[ilayer]) {
	303	ncontributors = ilayer2;
	304	pos[0] = xvtx;
	305	pos[1] = yvtx;
	306	pos[2] = zvtx;
	307	err[0]=TMath::Sqrt(0.25*(e2xclOutLay[i1InnLayBest]+e2xclOutLay[i2InnLayBest]));
	308	err[1]=TMath::Sqrt(0.25*(e2yclOutLay[i1InnLayBest]+e2yclOutLay[i2InnLayBest]));
	309	err[2]=TMath::Sqrt(0.25*(e2zclOutLay[i1InnLayBest]+e2zclOutLay[i2InnLayBest]));
	310	}
81ba01de	311	}
7e33c4e6	312	} // give it a try exchanging InnLay and OutLay
b107ee21	313
3acc14d5	314	fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
3acc14d5	315	fCurrentVertex->SetTitle("cosmics fake vertex");
b107ee21	316	fCurrentVertex->SetNContributors(ncontributors);
b107ee21	317	//fCurrentVertex->Print();
3acc14d5	318
	319	delete recpoints;
	320	itsLoader->UnloadRecPoints();
	321
	322	return fCurrentVertex;
	323	}
	324	//-------------------------------------------------------------------------
	325	void AliITSVertexerCosmics::FindVertices()
	326	{
	327	// computes the vertices of the events in the range FirstEvent - LastEvent
	328	AliRunLoader *rl = AliRunLoader::GetRunLoader();
	329	AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
	330	itsLoader->ReloadRecPoints();
	331	for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
	332	// printf("Processing event %d\n",i);
	333	rl->GetEvent(i);
	334	FindVertexForCurrentEvent(i);
	335	if(fCurrentVertex){
	336	WriteCurrentVertex();
	337	}
	338	}
	339	}
	340	//-------------------------------------------------------------------------
	341	void AliITSVertexerCosmics::PrintStatus() const
	342	{
	343	// Print current status
	344	printf("=======================================================\n");
b8ed1a92	345	printf(" fMaxDistOnOuterLayer: %f\n",fMaxDistOnOuterLayer);
b8ed1a92	346	printf(" fMaxVtxRadius[0]: %f\n",fMaxVtxRadius[0]);
3acc14d5	347	printf(" fMinDist2Vtxs: %f\n",fMinDist2Vtxs);
3acc14d5	348	printf("=======================================================\n");
	349	}
	350	//-------------------------------------------------------------------------