[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.
//
// -----------------------------------------------------------------
// 
// 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.
// 
//  
//  
//
//____________________________________________________________________

#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() {
  // 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];
  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.);
  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) {
  // Copy constructor
  // Copies are not allowed. The method is protected to avoid misuse.
  Error("AliITSMultReconstructor","Copy constructor not allowed\n");
}

//______________________________________________________________________
AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& /* mr */){
  // Assignment operator
  // Assignment is not allowed. The method is protected to avoid misuse.
  Error("= operator","Assignment operator not allowed\n");
  return *this;
}

//______________________________________________________________________
AliITSMultReconstructor::~AliITSMultReconstructor(){
  // Destructor
  if(fhClustersDPhiAcc)delete fhClustersDPhiAcc;
  if(fhClustersDThetaAcc)delete fhClustersDThetaAcc;
  if(fhClustersDZetaAcc)delete fhClustersDZetaAcc;
  if(fhClustersDPhiAll)delete fhClustersDPhiAll;
  if(fhClustersDThetaAll)delete fhClustersDThetaAll;
  if(fhClustersDZetaAll)delete fhClustersDZetaAll;
  if(fhDPhiVsDThetaAll)delete fhDPhiVsDThetaAll;
  if(fhDPhiVsDThetaAcc)delete fhDPhiVsDThetaAcc;
  if(fhDPhiVsDZetaAll)delete fhDPhiVsDZetaAll;
  if(fhDPhiVsDZetaAcc)delete fhDPhiVsDZetaAcc;
  if(fhetaTracklets)delete fhetaTracklets;
  if(fhphiTracklets)delete fhphiTracklets;
  if(fhetaClustersLay1)delete fhetaClustersLay1;
  if(fhphiClustersLay1)delete fhphiClustersLay1;
  if(fClustersLay1){
    for(Int_t i=0; i<300000; i++) {
      delete [] fClustersLay1[i];
      delete [] fClustersLay2[i];
      delete [] fTracklets[i];
    }
    delete fClustersLay1;
    delete  fClustersLay2;
    delete fTracklets;
  }
  if(fAssociationFlag)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; 

  // 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
       
      fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
      fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
      // Store a flag. This will indicate that the "tracklet" 
      // was indeed a single cluster! 
      fTracklets[fNTracklets][2] = -999999.;       
      AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)", 
		      fNTracklets, iC1));
      fNTracklets++;
    }

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