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

//____________________________________________________________________
// 
// AliITSMultReconstructor - find clusters in the pixels (theta and
// phi) and tracklets.
// 
// These can be used to extract charged particles multiplcicity from the ITS.
//
// A tracklet consist of two ITS clusters, one in the first pixel
// layer and one in the second. The clusters are associates if the 
// differencies in Phi (azimuth) and Zeta (longitudinal) are inside 
// a fiducial volume. In case of multiple candidates it is selected the
// candidate with minimum distance in Phi. 
// The parameter AssociationChoice allows to control if two clusters 
// in layer 2 can be associated to the same cluster in layer 1 or not.
//
// -----------------------------------------------------------------
// 
// TODO: 
// 
// - Introduce a rough pt estimation from the difference in phi ? 
// - Allow for a more refined selection criterium in case of multiple 
//   candidates (for instance by introducing weights for the difference 
//   in Phi and Zeta). 
//
//____________________________________________________________________

#include "AliITSMultReconstructor.h"

#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"


#include "AliITSRecPoint.h"
#include "AliITSgeom.h"
#include "AliLog.h"

//____________________________________________________________________
ClassImp(AliITSMultReconstructor)

//____________________________________________________________________
AliITSMultReconstructor::AliITSMultReconstructor() {

  fGeometry =0;

  SetHistOn();
  SetPhiWindow();
  SetZetaWindow();
  SetOnlyOneTrackletPerC2();

  fClustersLay1       = new Float_t*[300000];
  fClustersLay2       = new Float_t*[300000];
  fTracklets          = new Float_t*[300000];
  fAssociationFlag    = new Bool_t[300000];

  for(Int_t i=0; i<300000; i++) {
    fClustersLay1[i]       = new Float_t[3];
    fClustersLay2[i]       = new Float_t[3];
    fTracklets[i]          = new Float_t[3];
    fAssociationFlag[i]    = kFALSE;
  }

  // definition of histograms
  fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,-0.1,0.1);
  fhClustersDPhiAcc->SetDirectory(0);
  fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
  fhClustersDThetaAcc->SetDirectory(0);
  fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
  fhClustersDZetaAcc->SetDirectory(0);

  fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,-0.1,0.1);
  fhDPhiVsDZetaAcc->SetDirectory(0);
  fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
  fhDPhiVsDThetaAcc->SetDirectory(0);

  fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,-0.5,0.5);
  fhClustersDPhiAll->SetDirectory(0);
  fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
  fhClustersDThetaAll->SetDirectory(0);
  fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
  fhClustersDZetaAll->SetDirectory(0);

  fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,-0.5,0.5);
  fhDPhiVsDZetaAll->SetDirectory(0);
  fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,-0.5,0.5);
  fhDPhiVsDThetaAll->SetDirectory(0);

  fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
  fhetaTracklets->SetDirectory(0);
  fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100,-3.14159,3.14159);
  fhphiTracklets->SetDirectory(0);
  fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
  fhetaClustersLay1->SetDirectory(0);
  fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100,-3.141,3.141);
  fhphiClustersLay1->SetDirectory(0);
}
//____________________________________________________________________
AliITSMultReconstructor::~AliITSMultReconstructor() {
  // Destructor

  fGeometry = 0x0;
  for(Int_t i=0; i<300000; i++) {
    delete [] fClustersLay1[i];
    delete [] fClustersLay2[i];
    delete [] fTracklets[i];
  }
  
  delete [] fClustersLay1;
  delete [] fClustersLay2;
  delete [] fTracklets;
  delete [] fAssociationFlag;
  
  delete fhClustersDPhiAcc;
  delete fhClustersDThetaAcc;
  delete fhClustersDZetaAcc;
  delete fhClustersDPhiAll;
  delete fhClustersDThetaAll;
  delete fhClustersDZetaAll;
 
  delete fhDPhiVsDThetaAll;
  delete fhDPhiVsDThetaAcc;
  delete fhDPhiVsDZetaAll;
  delete fhDPhiVsDZetaAcc;

  delete fhetaTracklets;
  delete fhphiTracklets;
  delete fhetaClustersLay1;
  delete fhphiClustersLay1;
}

//____________________________________________________________________
void
AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
  //
  // - calls LoadClusterArray that finds the position of the clusters
  //   (in global coord) 
  // - convert the cluster coordinates to theta, phi (seen from the
  //   interaction vertex). The third coordinate is used for ....
  // - makes an array of tracklets 
  //   
  // After this method has been called, the clusters of the two layers
  // and the tracklets can be retrieved by calling the Get'er methods.

  // reset counters
  fNClustersLay1 = 0;
  fNClustersLay2 = 0;
  fNTracklets = 0; 

  // loading the clusters 
  LoadClusterArrays(clusterTree);
 
  // find the tracklets
  AliDebug(1,"Looking for tracklets... ");  

  //###########################################################
  // Loop on layer 1 : finding theta, phi and z 
  for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {    
    Float_t x = fClustersLay1[iC1][0] - vtx[0];
    Float_t y = fClustersLay1[iC1][1] - vtx[1];
    Float_t z = fClustersLay1[iC1][2] - vtx[2];

    Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
			       TMath::Power(y,2) +
			       TMath::Power(z,2));
    
    fClustersLay1[iC1][0] = TMath::ACos(z/r);  // Store Theta
    fClustersLay1[iC1][1] = TMath::ATan2(x,y);  // Store Phi
    fClustersLay1[iC1][2] = z/r;               // Store scaled z 
    if (fHistOn) {
      Float_t eta=fClustersLay1[iC1][0];
      eta= TMath::Tan(eta/2.);
      eta=-TMath::Log(eta);
      fhetaClustersLay1->Fill(eta);    
      fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);    
    }      
}
  
  // Loop on layer 2 : finding theta, phi and r   
  for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {    
    Float_t x = fClustersLay2[iC2][0] - vtx[0];
    Float_t y = fClustersLay2[iC2][1] - vtx[1];
    Float_t z = fClustersLay2[iC2][2] - vtx[2];
   
    Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
			       TMath::Power(y,2) +
			       TMath::Power(z,2));
    
    fClustersLay2[iC2][0] = TMath::ACos(z/r);  // Store Theta
    fClustersLay2[iC2][1] = TMath::ATan2(x,y);  // Store Phi
    fClustersLay2[iC2][2] = z;                 // Store z

 // this only needs to be initialized for the fNClustersLay2 first associations
    fAssociationFlag[iC2] = kFALSE;
  }  
  
  //###########################################################
  // Loop on layer 1 
  for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {    

    // reset of variables for multiple candidates
    Int_t  iC2WithBestDist = 0;     // reset 
    Float_t Distmin        = 100.;  // just to put a huge number! 
    Float_t dPhimin        = 0.;  // Used for histograms only! 
    Float_t dThetamin      = 0.;  // Used for histograms only! 
    Float_t dZetamin       = 0.;  // Used for histograms only! 
    
    // Loop on layer 2 
    for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {      
      
      // The following excludes double associations
      if (!fAssociationFlag[iC2]) {
	
	// find the difference in angles
	Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
	Float_t dPhi   = fClustersLay2[iC2][1] - fClustersLay1[iC1][1];
	
	// find the difference in z (between linear projection from layer 1
	// and the actual point: Dzeta= z1/r1*r2 -z2) 	
	Float_t r2   = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
        Float_t dZeta  = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2]; 

 	if (fHistOn) {
	  fhClustersDPhiAll->Fill(dPhi);    
	  fhClustersDThetaAll->Fill(dTheta);    
	  fhClustersDZetaAll->Fill(dZeta);    
	  fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
	  fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
	}
	// make "elliptical" cut in Phi and Zeta! 
	Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
	if (d>1) continue;      
	
	//look for the minimum distance: the minimum is in iC2WithBestDist
       	if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < Distmin ) {
	  Distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
	  dPhimin = dPhi;
	  dThetamin = dTheta;
	  dZetamin = dZeta; 
	  iC2WithBestDist = iC2;
	}
      } 
    } // end of loop over clusters in layer 2 
    
    if (Distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1

	if (fHistOn) {
	  fhClustersDPhiAcc->Fill(dPhimin);    
	  fhClustersDThetaAcc->Fill(dThetamin);    
	  fhClustersDZetaAcc->Fill(dZetamin);    
	  fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
	  fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
	}
      
      if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
      
      // store the tracklet
      
      // use the theta from the clusters in the first layer 
      fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
      // use the phi from the clusters in the first layer 
      fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
      // Store the difference between phi1 and phi2
      fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];       
  
 	if (fHistOn) {
	  Float_t eta=fTracklets[fNTracklets][0];
	  eta= TMath::Tan(eta/2.);
	  eta=-TMath::Log(eta);
	  fhetaTracklets->Fill(eta);    
	  fhphiTracklets->Fill(fTracklets[fNTracklets][1]);    
	}
      fNTracklets++;
      
      AliDebug(1,Form(" Adding tracklet candidate %d (cluster %d  of layer 1 and %d  of layer 2)", fNTracklets, iC1));
    }
  } // end of loop over clusters in layer 1
  
  AliDebug(1,Form("%d tracklets found", fNTracklets));
}

//____________________________________________________________________
void
AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
  // This method
  // - gets the clusters from the cluster tree 
  // - convert them into global coordinates 
  // - store them in the internal arrays
  
  AliDebug(1,"Loading clusters ...");
  
  fNClustersLay1 = 0;
  fNClustersLay2 = 0;
  
  TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
  TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");

  itsClusterBranch->SetAddress(&itsClusters);

  Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();  
 
  // loop over the its subdetectors
  for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
    
    if (!itsClusterTree->GetEvent(iIts)) 
      continue;
    
    Int_t nClusters = itsClusters->GetEntriesFast();
    
    // stuff needed to get the global coordinates
    Double_t rot[9];   fGeometry->GetRotMatrix(iIts,rot);
    Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det);
    Float_t tx,ty,tz;  fGeometry->GetTrans(lay,lad,det,tx,ty,tz);
    
    // Below:
    // "alpha" is the angle from the global X-axis to the
    //         local GEANT X'-axis  ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
    // "phi" is the angle from the global X-axis to the
    //       local cluster X"-axis
    
    Double_t alpha   = TMath::ATan2(rot[1],rot[0])+TMath::Pi();
    Double_t itsPhi = TMath::Pi()/2+alpha;
    
    if (lay==1) itsPhi+=TMath::Pi();
    Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi);
    Double_t r=tx*cp+ty*sp;
    
    // loop over clusters
    while(nClusters--) {
      AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);	
      
      if (cluster->GetLayer()>1) 
	continue;            
      
      Float_t x = r*cp - cluster->GetY()*sp;
      Float_t y = r*sp + cluster->GetY()*cp;
      Float_t z = cluster->GetZ();      
      
      if (cluster->GetLayer()==0) {
	fClustersLay1[fNClustersLay1][0] = x;
	fClustersLay1[fNClustersLay1][1] = y;
	fClustersLay1[fNClustersLay1][2] = z;
	fNClustersLay1++;
      }
      if (cluster->GetLayer()==1) {	
	fClustersLay2[fNClustersLay2][0] = x;
	fClustersLay2[fNClustersLay2][1] = y;
	fClustersLay2[fNClustersLay2][2] = z;
	fNClustersLay2++;
      }
      
    }// end of cluster loop
  } // end of its "subdetector" loop  
  
  AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
}
//____________________________________________________________________
void
AliITSMultReconstructor::SaveHists() {
  
  if (!fHistOn)
    return;

  fhClustersDPhiAll->Write();
  fhClustersDThetaAll->Write();
  fhClustersDZetaAll->Write();
  fhDPhiVsDThetaAll->Write();
  fhDPhiVsDZetaAll->Write();

  fhClustersDPhiAcc->Write();
  fhClustersDThetaAcc->Write();
  fhClustersDZetaAcc->Write();
  fhDPhiVsDThetaAcc->Write();
  fhDPhiVsDZetaAcc->Write();

  fhetaTracklets->Write();
  fhphiTracklets->Write();
  fhetaClustersLay1->Write();
  fhphiClustersLay1->Write();
}
Commit	Line	Data
ac903f1b	1	//____________________________________________________________________
	2	//
	3	// AliITSMultReconstructor - find clusters in the pixels (theta and
	4	// phi) and tracklets.
	5	//
	6	// These can be used to extract charged particles multiplcicity from the ITS.
	7	//
	8	// A tracklet consist of two ITS clusters, one in the first pixel
	9	// layer and one in the second. The clusters are associates if the
	10	// differencies in Phi (azimuth) and Zeta (longitudinal) are inside
	11	// a fiducial volume. In case of multiple candidates it is selected the
	12	// candidate with minimum distance in Phi.
	13	// The parameter AssociationChoice allows to control if two clusters
	14	// in layer 2 can be associated to the same cluster in layer 1 or not.
	15	//
	16	// -----------------------------------------------------------------
	17	//
	18	// TODO:
	19	//
	20	// - Introduce a rough pt estimation from the difference in phi ?
	21	// - Allow for a more refined selection criterium in case of multiple
	22	// candidates (for instance by introducing weights for the difference
	23	// in Phi and Zeta).
	24	//
	25	//____________________________________________________________________
	26
	27	#include "AliITSMultReconstructor.h"
	28
	29	#include "TTree.h"
	30	#include "TH1F.h"
	31	#include "TH2F.h"
	32
	33
b51872de	34	#include "AliITSRecPoint.h"
ac903f1b	35	#include "AliITSgeom.h"
	36	#include "AliLog.h"
	37
	38	//____________________________________________________________________
0762f3a8	39	ClassImp(AliITSMultReconstructor)
ac903f1b	40
	41	//____________________________________________________________________
	42	AliITSMultReconstructor::AliITSMultReconstructor() {
	43
	44	fGeometry =0;
	45
	46	SetHistOn();
	47	SetPhiWindow();
	48	SetZetaWindow();
	49	SetOnlyOneTrackletPerC2();
	50
	51	fClustersLay1 = new Float_t*[300000];
	52	fClustersLay2 = new Float_t*[300000];
	53	fTracklets = new Float_t*[300000];
	54	fAssociationFlag = new Bool_t[300000];
	55
	56	for(Int_t i=0; i<300000; i++) {
	57	fClustersLay1[i] = new Float_t[3];
	58	fClustersLay2[i] = new Float_t[3];
	59	fTracklets[i] = new Float_t[3];
	60	fAssociationFlag[i] = kFALSE;
	61	}
	62
	63	// definition of histograms
ddced3c8	64	fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,-0.1,0.1);
	65	fhClustersDPhiAcc->SetDirectory(0);
	66	fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
	67	fhClustersDThetaAcc->SetDirectory(0);
	68	fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
	69	fhClustersDZetaAcc->SetDirectory(0);
	70
	71	fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,-0.1,0.1);
	72	fhDPhiVsDZetaAcc->SetDirectory(0);
	73	fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
ac903f1b	74	fhDPhiVsDThetaAcc->SetDirectory(0);
ac903f1b	75
ddced3c8	76	fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,-0.5,0.5);
	77	fhClustersDPhiAll->SetDirectory(0);
	78	fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
	79	fhClustersDThetaAll->SetDirectory(0);
	80	fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
	81	fhClustersDZetaAll->SetDirectory(0);
	82
	83	fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,-0.5,0.5);
	84	fhDPhiVsDZetaAll->SetDirectory(0);
	85	fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,-0.5,0.5);
	86	fhDPhiVsDThetaAll->SetDirectory(0);
	87
	88	fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
	89	fhetaTracklets->SetDirectory(0);
	90	fhphiTracklets = new TH1F("phiTracklets", "phi", 100,-3.14159,3.14159);
	91	fhphiTracklets->SetDirectory(0);
	92	fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
	93	fhetaClustersLay1->SetDirectory(0);
	94	fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100,-3.141,3.141);
	95	fhphiClustersLay1->SetDirectory(0);
ac903f1b	96	}
ddced3c8	97	//____________________________________________________________________
	98	AliITSMultReconstructor::~AliITSMultReconstructor() {
	99	// Destructor
	100
	101	fGeometry = 0x0;
	102	for(Int_t i=0; i<300000; i++) {
	103	delete [] fClustersLay1[i];
	104	delete [] fClustersLay2[i];
	105	delete [] fTracklets[i];
	106	}
	107
	108	delete [] fClustersLay1;
	109	delete [] fClustersLay2;
	110	delete [] fTracklets;
	111	delete [] fAssociationFlag;
	112
	113	delete fhClustersDPhiAcc;
	114	delete fhClustersDThetaAcc;
	115	delete fhClustersDZetaAcc;
	116	delete fhClustersDPhiAll;
	117	delete fhClustersDThetaAll;
	118	delete fhClustersDZetaAll;
	119
	120	delete fhDPhiVsDThetaAll;
	121	delete fhDPhiVsDThetaAcc;
	122	delete fhDPhiVsDZetaAll;
	123	delete fhDPhiVsDZetaAcc;
ac903f1b	124
ddced3c8	125	delete fhetaTracklets;
	126	delete fhphiTracklets;
	127	delete fhetaClustersLay1;
	128	delete fhphiClustersLay1;
	129	}
ac903f1b	130
	131	//____________________________________________________________________
	132	void
	133	AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
	134	//
	135	// - calls LoadClusterArray that finds the position of the clusters
	136	// (in global coord)
	137	// - convert the cluster coordinates to theta, phi (seen from the
	138	// interaction vertex). The third coordinate is used for ....
	139	// - makes an array of tracklets
	140	//
	141	// After this method has been called, the clusters of the two layers
	142	// and the tracklets can be retrieved by calling the Get'er methods.
	143
ac903f1b	144	// reset counters
	145	fNClustersLay1 = 0;
	146	fNClustersLay2 = 0;
	147	fNTracklets = 0;
	148
	149	// loading the clusters
	150	LoadClusterArrays(clusterTree);
ddced3c8	151
ac903f1b	152	// find the tracklets
	153	AliDebug(1,"Looking for tracklets... ");
	154
	155	//###########################################################
	156	// Loop on layer 1 : finding theta, phi and z
	157	for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
	158	Float_t x = fClustersLay1[iC1][0] - vtx[0];
	159	Float_t y = fClustersLay1[iC1][1] - vtx[1];
	160	Float_t z = fClustersLay1[iC1][2] - vtx[2];
ddced3c8	161
ac903f1b	162	Float_t r = TMath::Sqrt(TMath::Power(x,2) +
	163	TMath::Power(y,2) +
	164	TMath::Power(z,2));
	165
	166	fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
ddced3c8	167	fClustersLay1[iC1][1] = TMath::ATan2(x,y); // Store Phi
ac903f1b	168	fClustersLay1[iC1][2] = z/r; // Store scaled z
ddced3c8	169	if (fHistOn) {
	170	Float_t eta=fClustersLay1[iC1][0];
	171	eta= TMath::Tan(eta/2.);
	172	eta=-TMath::Log(eta);
	173	fhetaClustersLay1->Fill(eta);
	174	fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
	175	}
	176	}
ac903f1b	177
	178	// Loop on layer 2 : finding theta, phi and r
	179	for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
	180	Float_t x = fClustersLay2[iC2][0] - vtx[0];
	181	Float_t y = fClustersLay2[iC2][1] - vtx[1];
	182	Float_t z = fClustersLay2[iC2][2] - vtx[2];
ddced3c8	183
ac903f1b	184	Float_t r = TMath::Sqrt(TMath::Power(x,2) +
	185	TMath::Power(y,2) +
	186	TMath::Power(z,2));
	187
	188	fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
ddced3c8	189	fClustersLay2[iC2][1] = TMath::ATan2(x,y); // Store Phi
ac903f1b	190	fClustersLay2[iC2][2] = z; // Store z
ac903f1b	191
ddced3c8	192	// this only needs to be initialized for the fNClustersLay2 first associations
ac903f1b	193	fAssociationFlag[iC2] = kFALSE;
	194	}
	195
	196	//###########################################################
	197	// Loop on layer 1
	198	for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
	199
	200	// reset of variables for multiple candidates
ddced3c8	201	Int_t iC2WithBestDist = 0; // reset
	202	Float_t Distmin = 100.; // just to put a huge number!
	203	Float_t dPhimin = 0.; // Used for histograms only!
	204	Float_t dThetamin = 0.; // Used for histograms only!
	205	Float_t dZetamin = 0.; // Used for histograms only!
ac903f1b	206
	207	// Loop on layer 2
	208	for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
	209
	210	// The following excludes double associations
	211	if (!fAssociationFlag[iC2]) {
	212
	213	// find the difference in angles
	214	Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
	215	Float_t dPhi = fClustersLay2[iC2][1] - fClustersLay1[iC1][1];
	216
	217	// find the difference in z (between linear projection from layer 1
	218	// and the actual point: Dzeta= z1/r1*r2 -z2)
ddced3c8	219	Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
	220	Float_t dZeta = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
	221
	222	if (fHistOn) {
	223	fhClustersDPhiAll->Fill(dPhi);
	224	fhClustersDThetaAll->Fill(dTheta);
	225	fhClustersDZetaAll->Fill(dZeta);
ac903f1b	226	fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
ddced3c8	227	fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
ac903f1b	228	}
	229	// make "elliptical" cut in Phi and Zeta!
	230	Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
	231	if (d>1) continue;
	232
ddced3c8	233	//look for the minimum distance: the minimum is in iC2WithBestDist
	234	if (TMath::Sqrt(dZetadZeta+(r2dPhir2dPhi)) < Distmin ) {
	235	Distmin=TMath::Sqrt(dZetadZeta + (r2dPhir2dPhi));
	236	dPhimin = dPhi;
	237	dThetamin = dTheta;
	238	dZetamin = dZeta;
	239	iC2WithBestDist = iC2;
ac903f1b	240	}
	241	}
	242	} // end of loop over clusters in layer 2
	243
ddced3c8	244	if (Distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
	245
	246	if (fHistOn) {
	247	fhClustersDPhiAcc->Fill(dPhimin);
	248	fhClustersDThetaAcc->Fill(dThetamin);
	249	fhClustersDZetaAcc->Fill(dZetamin);
	250	fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
	251	fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
	252	}
ac903f1b	253
ddced3c8	254	if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
ac903f1b	255
	256	// store the tracklet
	257
ddced3c8	258	// use the theta from the clusters in the first layer
ddced3c8	259	fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
ac903f1b	260	// use the phi from the clusters in the first layer
	261	fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
	262	// Store the difference between phi1 and phi2
ddced3c8	263	fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
	264
	265	if (fHistOn) {
	266	Float_t eta=fTracklets[fNTracklets][0];
	267	eta= TMath::Tan(eta/2.);
	268	eta=-TMath::Log(eta);
	269	fhetaTracklets->Fill(eta);
	270	fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
	271	}
ac903f1b	272	fNTracklets++;
	273
	274	AliDebug(1,Form(" Adding tracklet candidate %d (cluster %d of layer 1 and %d of layer 2)", fNTracklets, iC1));
	275	}
	276	} // end of loop over clusters in layer 1
	277
	278	AliDebug(1,Form("%d tracklets found", fNTracklets));
	279	}
	280
	281	//____________________________________________________________________
	282	void
	283	AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
	284	// This method
	285	// - gets the clusters from the cluster tree
	286	// - convert them into global coordinates
	287	// - store them in the internal arrays
	288
	289	AliDebug(1,"Loading clusters ...");
	290
	291	fNClustersLay1 = 0;
	292	fNClustersLay2 = 0;
	293
b51872de	294	TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
b51872de	295	TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
ddced3c8	296
ac903f1b	297	itsClusterBranch->SetAddress(&itsClusters);
ddced3c8	298
ac903f1b	299	Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
ddced3c8	300
ac903f1b	301	// loop over the its subdetectors
	302	for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
	303
	304	if (!itsClusterTree->GetEvent(iIts))
	305	continue;
	306
	307	Int_t nClusters = itsClusters->GetEntriesFast();
	308
	309	// stuff needed to get the global coordinates
	310	Double_t rot[9]; fGeometry->GetRotMatrix(iIts,rot);
	311	Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det);
	312	Float_t tx,ty,tz; fGeometry->GetTrans(lay,lad,det,tx,ty,tz);
	313
	314	// Below:
	315	// "alpha" is the angle from the global X-axis to the
	316	// local GEANT X'-axis ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
	317	// "phi" is the angle from the global X-axis to the
	318	// local cluster X"-axis
	319
	320	Double_t alpha = TMath::ATan2(rot[1],rot[0])+TMath::Pi();
	321	Double_t itsPhi = TMath::Pi()/2+alpha;
	322
	323	if (lay==1) itsPhi+=TMath::Pi();
	324	Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi);
	325	Double_t r=txcp+tysp;
	326
	327	// loop over clusters
	328	while(nClusters--) {
b51872de	329	AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
ac903f1b	330
	331	if (cluster->GetLayer()>1)
	332	continue;
	333
	334	Float_t x = rcp - cluster->GetY()sp;
	335	Float_t y = rsp + cluster->GetY()cp;
	336	Float_t z = cluster->GetZ();
	337
	338	if (cluster->GetLayer()==0) {
	339	fClustersLay1[fNClustersLay1][0] = x;
	340	fClustersLay1[fNClustersLay1][1] = y;
	341	fClustersLay1[fNClustersLay1][2] = z;
	342	fNClustersLay1++;
	343	}
	344	if (cluster->GetLayer()==1) {
	345	fClustersLay2[fNClustersLay2][0] = x;
	346	fClustersLay2[fNClustersLay2][1] = y;
	347	fClustersLay2[fNClustersLay2][2] = z;
	348	fNClustersLay2++;
	349	}
	350
	351	}// end of cluster loop
	352	} // end of its "subdetector" loop
	353
	354	AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
	355	}
	356	//____________________________________________________________________
	357	void
	358	AliITSMultReconstructor::SaveHists() {
	359
	360	if (!fHistOn)
	361	return;
	362
ddced3c8	363	fhClustersDPhiAll->Write();
	364	fhClustersDThetaAll->Write();
	365	fhClustersDZetaAll->Write();
ac903f1b	366	fhDPhiVsDThetaAll->Write();
ddced3c8	367	fhDPhiVsDZetaAll->Write();
	368
	369	fhClustersDPhiAcc->Write();
	370	fhClustersDThetaAcc->Write();
	371	fhClustersDZetaAcc->Write();
ac903f1b	372	fhDPhiVsDThetaAcc->Write();
ddced3c8	373	fhDPhiVsDZetaAcc->Write();
	374
	375	fhetaTracklets->Write();
	376	fhphiTracklets->Write();
	377	fhetaClustersLay1->Write();
	378	fhphiClustersLay1->Write();
ac903f1b	379	}