[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():
fGeometry(0),
fClustersLay1(0),
fClustersLay2(0),
fTracklets(0),
fAssociationFlag(0),
fNClustersLay1(0),
fNClustersLay2(0),
fNTracklets(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];
  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),
fGeometry(mr.fGeometry),
fClustersLay1(mr.fClustersLay1),
fClustersLay2(mr.fClustersLay2),
fTracklets(mr.fTracklets),
fAssociationFlag(mr.fAssociationFlag),
fNClustersLay1(mr.fNClustersLay1),
fNClustersLay2(mr.fNClustersLay2),
fNTracklets(mr.fNTracklets),
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 [] fClustersLay1;
  delete [] fClustersLay2;
  delete [] fTracklets;

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