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

/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/

/* $Id$ */

//____________________________________________________________________
// 
// 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.
//
// Two methods return the number of traklets and the number of clusters 
// in the first SPD layer (GetNTracklets GetNSingleClusters)
//
// -----------------------------------------------------------------
// 
// NOTE: The cuts on phi and zeta depends on the interacting system (p-p  
//  or Pb-Pb). Please, check the file AliITSMultReconstructor.h and be 
//  sure that SetPhiWindow and SetZetaWindow are defined accordingly.
// 
//  Author :  Tiziano Virgili 
//  
//  
//
//____________________________________________________________________

#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TTree.h>

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

//____________________________________________________________________
ClassImp(AliITSMultReconstructor)


//____________________________________________________________________
AliITSMultReconstructor::AliITSMultReconstructor():
fGeometry(0),
fClustersLay1(0),
fClustersLay2(0),
fTracklets(0),
fSClusters(0),
fAssociationFlag(0),
fNClustersLay1(0),
fNClustersLay2(0),
fNTracklets(0),
fNSingleCluster(0),
fPhiWindow(0),
fZetaWindow(0),
fOnlyOneTrackletPerC2(0),
fHistOn(0),
fhClustersDPhiAcc(0),
fhClustersDThetaAcc(0),
fhClustersDZetaAcc(0),
fhClustersDPhiAll(0),
fhClustersDThetaAll(0),
fhClustersDZetaAll(0),
fhDPhiVsDThetaAll(0),
fhDPhiVsDThetaAcc(0),
fhDPhiVsDZetaAll(0),
fhDPhiVsDZetaAcc(0),
fhetaTracklets(0),
fhphiTracklets(0),
fhetaClustersLay1(0),
fhphiClustersLay1(0){
  // Method to reconstruct the charged particles multiplicity with the 
  // SPD (tracklets).

  fGeometry =0;

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

  fClustersLay1       = new Float_t*[300000];
  fClustersLay2       = new Float_t*[300000];
  fTracklets          = new Float_t*[300000];
  fSClusters          = 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];
    fSClusters[i]           = new Float_t[2];
    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.);
  fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100,-3.14159,3.14159);
  fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
  fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100,-3.141,3.141);

}

//______________________________________________________________________
AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
fGeometry(mr.fGeometry),
fClustersLay1(mr.fClustersLay1),
fClustersLay2(mr.fClustersLay2),
fTracklets(mr.fTracklets),
fSClusters(mr.fSClusters),
fAssociationFlag(mr.fAssociationFlag),
fNClustersLay1(mr.fNClustersLay1),
fNClustersLay2(mr.fNClustersLay2),
fNTracklets(mr.fNTracklets),
fNSingleCluster(mr.fNSingleCluster),
fPhiWindow(mr.fPhiWindow),
fZetaWindow(mr.fZetaWindow),
fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
fHistOn(mr.fHistOn),
fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
fhClustersDPhiAll(mr.fhClustersDPhiAll),
fhClustersDThetaAll(mr.fhClustersDThetaAll),
fhClustersDZetaAll(mr.fhClustersDZetaAll),
fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
fhetaTracklets(mr.fhetaTracklets),
fhphiTracklets(mr.fhphiTracklets),
fhetaClustersLay1(mr.fhetaClustersLay1),
fhphiClustersLay1(mr.fhphiClustersLay1) {
  // Copy constructor

}

//______________________________________________________________________
AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
  // Assignment operator
  this->~AliITSMultReconstructor();
  new(this) AliITSMultReconstructor(mr);
  return *this;
}

//______________________________________________________________________
AliITSMultReconstructor::~AliITSMultReconstructor(){
  // Destructor

  // delete histograms
  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;

  // delete arrays
  for(Int_t i=0; i<300000; i++) {
    delete [] fClustersLay1[i];
    delete [] fClustersLay2[i];
    delete [] fTracklets[i];
    delete [] fSClusters[i];
  }
  delete [] fClustersLay1;
  delete [] fClustersLay2;
  delete [] fTracklets;
  delete [] fSClusters;

  delete [] fAssociationFlag;
}

//____________________________________________________________________
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; 
  fNSingleCluster = 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]);    
      }
      
      AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
      AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1, 
		      iC2WithBestDist));
      fNTracklets++;
    }

    // Delete the following else if you do not want to save Clusters! 

    else { // This means that the cluster has not been associated 

      // store the cluster
       
      fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
      fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
      AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)", 
		      fNSingleCluster, iC1));
      fNSingleCluster++;
    }

  } // 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() {
  // This method save the histograms on the output file
  // (only if fHistOn is TRUE). 
  
  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
7ca4655f	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
ac903f1b	18	//____________________________________________________________________
	19	//
	20	// AliITSMultReconstructor - find clusters in the pixels (theta and
	21	// phi) and tracklets.
	22	//
	23	// These can be used to extract charged particles multiplcicity from the ITS.
	24	//
	25	// A tracklet consist of two ITS clusters, one in the first pixel
	26	// layer and one in the second. The clusters are associates if the
	27	// differencies in Phi (azimuth) and Zeta (longitudinal) are inside
	28	// a fiducial volume. In case of multiple candidates it is selected the
	29	// candidate with minimum distance in Phi.
	30	// The parameter AssociationChoice allows to control if two clusters
	31	// in layer 2 can be associated to the same cluster in layer 1 or not.
	32	//
968e8539	33	// Two methods return the number of traklets and the number of clusters
	34	// in the first SPD layer (GetNTracklets GetNSingleClusters)
	35	//
ac903f1b	36	// -----------------------------------------------------------------
ac903f1b	37	//
3ef75756	38	// NOTE: The cuts on phi and zeta depends on the interacting system (p-p
	39	// or Pb-Pb). Please, check the file AliITSMultReconstructor.h and be
	40	// sure that SetPhiWindow and SetZetaWindow are defined accordingly.
ac903f1b	41	//
968e8539	42	// Author : Tiziano Virgili
3ef75756	43	//
3ef75756	44	//
ac903f1b	45	//
	46	//____________________________________________________________________
	47
7ca4655f	48	#include <TClonesArray.h>
	49	#include <TH1F.h>
	50	#include <TH2F.h>
	51	#include <TTree.h>
ac903f1b	52
7ca4655f	53	#include "AliITSMultReconstructor.h"
b51872de	54	#include "AliITSRecPoint.h"
ac903f1b	55	#include "AliITSgeom.h"
	56	#include "AliLog.h"
	57
	58	//____________________________________________________________________
0762f3a8	59	ClassImp(AliITSMultReconstructor)
ac903f1b	60
3ef75756	61
ac903f1b	62	//____________________________________________________________________
7537d03c	63	AliITSMultReconstructor::AliITSMultReconstructor():
	64	fGeometry(0),
	65	fClustersLay1(0),
	66	fClustersLay2(0),
	67	fTracklets(0),
968e8539	68	fSClusters(0),
7537d03c	69	fAssociationFlag(0),
	70	fNClustersLay1(0),
	71	fNClustersLay2(0),
	72	fNTracklets(0),
968e8539	73	fNSingleCluster(0),
7537d03c	74	fPhiWindow(0),
	75	fZetaWindow(0),
	76	fOnlyOneTrackletPerC2(0),
	77	fHistOn(0),
	78	fhClustersDPhiAcc(0),
	79	fhClustersDThetaAcc(0),
	80	fhClustersDZetaAcc(0),
	81	fhClustersDPhiAll(0),
	82	fhClustersDThetaAll(0),
	83	fhClustersDZetaAll(0),
	84	fhDPhiVsDThetaAll(0),
	85	fhDPhiVsDThetaAcc(0),
	86	fhDPhiVsDZetaAll(0),
	87	fhDPhiVsDZetaAcc(0),
	88	fhetaTracklets(0),
	89	fhphiTracklets(0),
	90	fhetaClustersLay1(0),
	91	fhphiClustersLay1(0){
3ef75756	92	// Method to reconstruct the charged particles multiplicity with the
3ef75756	93	// SPD (tracklets).
ac903f1b	94
	95	fGeometry =0;
	96
	97	SetHistOn();
	98	SetPhiWindow();
	99	SetZetaWindow();
	100	SetOnlyOneTrackletPerC2();
	101
	102	fClustersLay1 = new Float_t*[300000];
	103	fClustersLay2 = new Float_t*[300000];
	104	fTracklets = new Float_t*[300000];
968e8539	105	fSClusters = new Float_t*[300000];
ac903f1b	106	fAssociationFlag = new Bool_t[300000];
	107
	108	for(Int_t i=0; i<300000; i++) {
	109	fClustersLay1[i] = new Float_t[3];
	110	fClustersLay2[i] = new Float_t[3];
	111	fTracklets[i] = new Float_t[3];
968e8539	112	fSClusters[i] = new Float_t[2];
ac903f1b	113	fAssociationFlag[i] = kFALSE;
	114	}
	115
	116	// definition of histograms
ddced3c8	117	fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,-0.1,0.1);
	118	fhClustersDPhiAcc->SetDirectory(0);
	119	fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
	120	fhClustersDThetaAcc->SetDirectory(0);
	121	fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
	122	fhClustersDZetaAcc->SetDirectory(0);
	123
	124	fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,-0.1,0.1);
	125	fhDPhiVsDZetaAcc->SetDirectory(0);
	126	fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
ac903f1b	127	fhDPhiVsDThetaAcc->SetDirectory(0);
ac903f1b	128
ddced3c8	129	fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,-0.5,0.5);
	130	fhClustersDPhiAll->SetDirectory(0);
	131	fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
	132	fhClustersDThetaAll->SetDirectory(0);
	133	fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
	134	fhClustersDZetaAll->SetDirectory(0);
	135
	136	fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,-0.5,0.5);
	137	fhDPhiVsDZetaAll->SetDirectory(0);
	138	fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,-0.5,0.5);
	139	fhDPhiVsDThetaAll->SetDirectory(0);
	140
	141	fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
ddced3c8	142	fhphiTracklets = new TH1F("phiTracklets", "phi", 100,-3.14159,3.14159);
ddced3c8	143	fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
ddced3c8	144	fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100,-3.141,3.141);
3ef75756	145
ac903f1b	146	}
ddced3c8	147
3ef75756	148	//______________________________________________________________________
7537d03c	149	AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
	150	fGeometry(mr.fGeometry),
	151	fClustersLay1(mr.fClustersLay1),
	152	fClustersLay2(mr.fClustersLay2),
	153	fTracklets(mr.fTracklets),
968e8539	154	fSClusters(mr.fSClusters),
7537d03c	155	fAssociationFlag(mr.fAssociationFlag),
	156	fNClustersLay1(mr.fNClustersLay1),
	157	fNClustersLay2(mr.fNClustersLay2),
	158	fNTracklets(mr.fNTracklets),
968e8539	159	fNSingleCluster(mr.fNSingleCluster),
7537d03c	160	fPhiWindow(mr.fPhiWindow),
	161	fZetaWindow(mr.fZetaWindow),
	162	fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
	163	fHistOn(mr.fHistOn),
	164	fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
	165	fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
	166	fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
	167	fhClustersDPhiAll(mr.fhClustersDPhiAll),
	168	fhClustersDThetaAll(mr.fhClustersDThetaAll),
	169	fhClustersDZetaAll(mr.fhClustersDZetaAll),
	170	fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
	171	fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
	172	fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
	173	fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
	174	fhetaTracklets(mr.fhetaTracklets),
	175	fhphiTracklets(mr.fhphiTracklets),
	176	fhetaClustersLay1(mr.fhetaClustersLay1),
	177	fhphiClustersLay1(mr.fhphiClustersLay1) {
3ef75756	178	// Copy constructor
7537d03c	179
3ef75756	180	}
	181
	182	//______________________________________________________________________
7537d03c	183	AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
3ef75756	184	// Assignment operator
7537d03c	185	this->~AliITSMultReconstructor();
7537d03c	186	new(this) AliITSMultReconstructor(mr);
3ef75756	187	return *this;
	188	}
	189
	190	//______________________________________________________________________
	191	AliITSMultReconstructor::~AliITSMultReconstructor(){
	192	// Destructor
1ba5b31c	193
	194	// delete histograms
	195	delete fhClustersDPhiAcc;
	196	delete fhClustersDThetaAcc;
	197	delete fhClustersDZetaAcc;
	198	delete fhClustersDPhiAll;
	199	delete fhClustersDThetaAll;
	200	delete fhClustersDZetaAll;
	201	delete fhDPhiVsDThetaAll;
	202	delete fhDPhiVsDThetaAcc;
	203	delete fhDPhiVsDZetaAll;
	204	delete fhDPhiVsDZetaAcc;
	205	delete fhetaTracklets;
	206	delete fhphiTracklets;
	207	delete fhetaClustersLay1;
	208	delete fhphiClustersLay1;
	209
	210	// delete arrays
	211	for(Int_t i=0; i<300000; i++) {
	212	delete [] fClustersLay1[i];
	213	delete [] fClustersLay2[i];
	214	delete [] fTracklets[i];
968e8539	215	delete [] fSClusters[i];
ddced3c8	216	}
1ba5b31c	217	delete [] fClustersLay1;
	218	delete [] fClustersLay2;
	219	delete [] fTracklets;
968e8539	220	delete [] fSClusters;
1ba5b31c	221
1ba5b31c	222	delete [] fAssociationFlag;
ddced3c8	223	}
ac903f1b	224
	225	//____________________________________________________________________
	226	void
	227	AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
	228	//
	229	// - calls LoadClusterArray that finds the position of the clusters
	230	// (in global coord)
	231	// - convert the cluster coordinates to theta, phi (seen from the
	232	// interaction vertex). The third coordinate is used for ....
	233	// - makes an array of tracklets
	234	//
	235	// After this method has been called, the clusters of the two layers
	236	// and the tracklets can be retrieved by calling the Get'er methods.
	237
ac903f1b	238	// reset counters
	239	fNClustersLay1 = 0;
	240	fNClustersLay2 = 0;
	241	fNTracklets = 0;
968e8539	242	fNSingleCluster = 0;
ac903f1b	243	// loading the clusters
ac903f1b	244	LoadClusterArrays(clusterTree);
3ef75756	245
ac903f1b	246	// find the tracklets
	247	AliDebug(1,"Looking for tracklets... ");
	248
	249	//###########################################################
	250	// Loop on layer 1 : finding theta, phi and z
	251	for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
	252	Float_t x = fClustersLay1[iC1][0] - vtx[0];
	253	Float_t y = fClustersLay1[iC1][1] - vtx[1];
	254	Float_t z = fClustersLay1[iC1][2] - vtx[2];
ddced3c8	255
ac903f1b	256	Float_t r = TMath::Sqrt(TMath::Power(x,2) +
	257	TMath::Power(y,2) +
	258	TMath::Power(z,2));
	259
	260	fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
ddced3c8	261	fClustersLay1[iC1][1] = TMath::ATan2(x,y); // Store Phi
ac903f1b	262	fClustersLay1[iC1][2] = z/r; // Store scaled z
ddced3c8	263	if (fHistOn) {
	264	Float_t eta=fClustersLay1[iC1][0];
	265	eta= TMath::Tan(eta/2.);
	266	eta=-TMath::Log(eta);
	267	fhetaClustersLay1->Fill(eta);
	268	fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
	269	}
	270	}
ac903f1b	271
	272	// Loop on layer 2 : finding theta, phi and r
	273	for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
	274	Float_t x = fClustersLay2[iC2][0] - vtx[0];
	275	Float_t y = fClustersLay2[iC2][1] - vtx[1];
	276	Float_t z = fClustersLay2[iC2][2] - vtx[2];
ddced3c8	277
ac903f1b	278	Float_t r = TMath::Sqrt(TMath::Power(x,2) +
	279	TMath::Power(y,2) +
	280	TMath::Power(z,2));
	281
	282	fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
ddced3c8	283	fClustersLay2[iC2][1] = TMath::ATan2(x,y); // Store Phi
ac903f1b	284	fClustersLay2[iC2][2] = z; // Store z
ac903f1b	285
ddced3c8	286	// this only needs to be initialized for the fNClustersLay2 first associations
ac903f1b	287	fAssociationFlag[iC2] = kFALSE;
	288	}
	289
	290	//###########################################################
	291	// Loop on layer 1
	292	for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
	293
	294	// reset of variables for multiple candidates
ddced3c8	295	Int_t iC2WithBestDist = 0; // reset
3ef75756	296	Float_t distmin = 100.; // just to put a huge number!
ddced3c8	297	Float_t dPhimin = 0.; // Used for histograms only!
	298	Float_t dThetamin = 0.; // Used for histograms only!
	299	Float_t dZetamin = 0.; // Used for histograms only!
ac903f1b	300
	301	// Loop on layer 2
	302	for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
	303
	304	// The following excludes double associations
	305	if (!fAssociationFlag[iC2]) {
	306
	307	// find the difference in angles
	308	Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
	309	Float_t dPhi = fClustersLay2[iC2][1] - fClustersLay1[iC1][1];
	310
	311	// find the difference in z (between linear projection from layer 1
	312	// and the actual point: Dzeta= z1/r1*r2 -z2)
ddced3c8	313	Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
	314	Float_t dZeta = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
	315
	316	if (fHistOn) {
	317	fhClustersDPhiAll->Fill(dPhi);
	318	fhClustersDThetaAll->Fill(dTheta);
	319	fhClustersDZetaAll->Fill(dZeta);
ac903f1b	320	fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
ddced3c8	321	fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
ac903f1b	322	}
	323	// make "elliptical" cut in Phi and Zeta!
	324	Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
3ef75756	325
ac903f1b	326	if (d>1) continue;
ac903f1b	327
ddced3c8	328	//look for the minimum distance: the minimum is in iC2WithBestDist
3ef75756	329	if (TMath::Sqrt(dZetadZeta+(r2dPhir2dPhi)) < distmin ) {
3ef75756	330	distmin=TMath::Sqrt(dZetadZeta + (r2dPhir2dPhi));
ddced3c8	331	dPhimin = dPhi;
	332	dThetamin = dTheta;
	333	dZetamin = dZeta;
	334	iC2WithBestDist = iC2;
ac903f1b	335	}
	336	}
	337	} // end of loop over clusters in layer 2
	338
3ef75756	339	if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
	340
	341	if (fHistOn) {
	342	fhClustersDPhiAcc->Fill(dPhimin);
	343	fhClustersDThetaAcc->Fill(dThetamin);
	344	fhClustersDZetaAcc->Fill(dZetamin);
	345	fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
	346	fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
	347	}
ac903f1b	348
ddced3c8	349	if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
ac903f1b	350
	351	// store the tracklet
	352
ddced3c8	353	// use the theta from the clusters in the first layer
ddced3c8	354	fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
ac903f1b	355	// use the phi from the clusters in the first layer
	356	fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
	357	// Store the difference between phi1 and phi2
ddced3c8	358	fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
ddced3c8	359
3ef75756	360	if (fHistOn) {
	361	Float_t eta=fTracklets[fNTracklets][0];
	362	eta= TMath::Tan(eta/2.);
	363	eta=-TMath::Log(eta);
	364	fhetaTracklets->Fill(eta);
	365	fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
	366	}
ac903f1b	367
3ef75756	368	AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
	369	AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1,
	370	iC2WithBestDist));
	371	fNTracklets++;
ac903f1b	372	}
3ef75756	373
	374	// Delete the following else if you do not want to save Clusters!
	375
	376	else { // This means that the cluster has not been associated
	377
	378	// store the cluster
	379
968e8539	380	fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
968e8539	381	fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
3ef75756	382	AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
968e8539	383	fNSingleCluster, iC1));
968e8539	384	fNSingleCluster++;
3ef75756	385	}
3ef75756	386
ac903f1b	387	} // end of loop over clusters in layer 1
	388
	389	AliDebug(1,Form("%d tracklets found", fNTracklets));
	390	}
	391
	392	//____________________________________________________________________
	393	void
	394	AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
	395	// This method
	396	// - gets the clusters from the cluster tree
	397	// - convert them into global coordinates
	398	// - store them in the internal arrays
	399
	400	AliDebug(1,"Loading clusters ...");
	401
	402	fNClustersLay1 = 0;
	403	fNClustersLay2 = 0;
	404
b51872de	405	TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
b51872de	406	TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
ddced3c8	407
ac903f1b	408	itsClusterBranch->SetAddress(&itsClusters);
ddced3c8	409
ac903f1b	410	Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
ddced3c8	411
ac903f1b	412	// loop over the its subdetectors
	413	for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
	414
	415	if (!itsClusterTree->GetEvent(iIts))
	416	continue;
	417
	418	Int_t nClusters = itsClusters->GetEntriesFast();
	419
	420	// stuff needed to get the global coordinates
	421	Double_t rot[9]; fGeometry->GetRotMatrix(iIts,rot);
	422	Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det);
	423	Float_t tx,ty,tz; fGeometry->GetTrans(lay,lad,det,tx,ty,tz);
	424
	425	// Below:
	426	// "alpha" is the angle from the global X-axis to the
	427	// local GEANT X'-axis ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
	428	// "phi" is the angle from the global X-axis to the
	429	// local cluster X"-axis
	430
	431	Double_t alpha = TMath::ATan2(rot[1],rot[0])+TMath::Pi();
	432	Double_t itsPhi = TMath::Pi()/2+alpha;
	433
	434	if (lay==1) itsPhi+=TMath::Pi();
	435	Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi);
	436	Double_t r=txcp+tysp;
	437
	438	// loop over clusters
	439	while(nClusters--) {
b51872de	440	AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
ac903f1b	441
	442	if (cluster->GetLayer()>1)
	443	continue;
	444
	445	Float_t x = rcp - cluster->GetY()sp;
	446	Float_t y = rsp + cluster->GetY()cp;
	447	Float_t z = cluster->GetZ();
	448
	449	if (cluster->GetLayer()==0) {
	450	fClustersLay1[fNClustersLay1][0] = x;
	451	fClustersLay1[fNClustersLay1][1] = y;
	452	fClustersLay1[fNClustersLay1][2] = z;
	453	fNClustersLay1++;
	454	}
	455	if (cluster->GetLayer()==1) {
	456	fClustersLay2[fNClustersLay2][0] = x;
	457	fClustersLay2[fNClustersLay2][1] = y;
	458	fClustersLay2[fNClustersLay2][2] = z;
	459	fNClustersLay2++;
	460	}
	461
	462	}// end of cluster loop
	463	} // end of its "subdetector" loop
	464
	465	AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
	466	}
	467	//____________________________________________________________________
	468	void
	469	AliITSMultReconstructor::SaveHists() {
3ef75756	470	// This method save the histograms on the output file
3ef75756	471	// (only if fHistOn is TRUE).
ac903f1b	472
	473	if (!fHistOn)
	474	return;
	475
ddced3c8	476	fhClustersDPhiAll->Write();
	477	fhClustersDThetaAll->Write();
	478	fhClustersDZetaAll->Write();
ac903f1b	479	fhDPhiVsDThetaAll->Write();
ddced3c8	480	fhDPhiVsDZetaAll->Write();
	481
	482	fhClustersDPhiAcc->Write();
	483	fhClustersDThetaAcc->Write();
	484	fhClustersDZetaAcc->Write();
ac903f1b	485	fhDPhiVsDThetaAcc->Write();
ddced3c8	486	fhDPhiVsDZetaAcc->Write();
	487
	488	fhetaTracklets->Write();
	489	fhphiTracklets->Write();
	490	fhetaClustersLay1->Write();
	491	fhphiClustersLay1->Write();
ac903f1b	492	}