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

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

#include "AliITSMultReconstructor.h"

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


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

//____________________________________________________________________
ClassImp(AliITSMultReconstructor)

//____________________________________________________________________
AliITSMultReconstructor::AliITSMultReconstructor() {

  fGeometry =0;

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

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

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

  // definition of histograms
  fhClustersDPhi   = new TH1F("dphi",  "dphi",  200,-0.1,0.1);
  fhClustersDPhi->SetDirectory(0);
  fhClustersDTheta = new TH1F("dtheta","dtheta",200,-0.1,0.1);
  fhClustersDTheta->SetDirectory(0);
  fhClustersDZeta = new TH1F("dzeta","dzeta",200,-0.2,0.2);
  fhClustersDZeta->SetDirectory(0);

  fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",200,-0.1,0.1,200,-0.1,0.1);
  fhDPhiVsDThetaAll->SetDirectory(0);
  fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",200,-0.1,0.1,200,-0.1,0.1);
  fhDPhiVsDThetaAcc->SetDirectory(0);

}


//____________________________________________________________________
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::ATan(y/x);  // Store Phi
    fClustersLay1[iC1][2] = z/r;               // Store scaled z 
  }
  
  // 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::ATan(y/x);  // 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   iC2WithBestPhi = 0;     // reset 
    Float_t dPhimin        = 100.;  // just to put a huge number! 
    
    // 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  = fClustersLay2[iC1][2]*r2 - fClustersLay2[iC2][2]; 
	
	if (fHistOn) {
	  fhClustersDPhi->Fill(dPhi);    
	  fhClustersDTheta->Fill(dTheta);    
	  fhClustersDZeta->Fill(dZeta);    
	  fhDPhiVsDThetaAll->Fill(dTheta, 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 in Phi: the minimum is in iC2WithBestPhi
	if (TMath::Abs(dPhi) < dPhimin) {
	  dPhimin = TMath::Abs(dPhi);
	  iC2WithBestPhi = iC2;
	}
      } 
    } // end of loop over clusters in layer 2 
    
    if (dPhimin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
      
      if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestPhi] = kTRUE; // flag the association
      
      // store the tracklet
      
      // use the average theta from the clusters in the two layers
      fTracklets[fNTracklets][0] = 0.5*(fClustersLay1[iC1][0]+fClustersLay2[iC2WithBestPhi][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[iC2WithBestPhi][1];         
      fNTracklets++;
      
      AliDebug(1,Form(" Adding tracklet candidate %d (cluster %d  of layer 1 and %d  of layer 2)", fNTracklets, iC1));
    }
  } // end of loop over clusters in layer 1
  
  AliDebug(1,Form("%d tracklets found", fNTracklets));
}

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

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