[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 "AliITSclusterV2.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.

  cout << " HEEEEEEEEEEEEEEEEEEEEE " << flush << endl;

  // 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("AliITSclusterV2");
  TBranch* itsClusterBranch=itsClusterTree->GetBranch("Clusters");
  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--) {
      AliITSclusterV2* cluster = (AliITSclusterV2*)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;

  cout << "Saving histograms" << endl;

  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
	34	#include "AliITSclusterV2.h"
	35	#include "AliITSgeom.h"
	36	#include "AliLog.h"
	37
	38	//____________________________________________________________________
	39	ClassImp(AliITSMultReconstructor);
	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
92	cout << " HEEEEEEEEEEEEEEEEEEEEE " << flush << endl;
93
94	// reset counters
95	fNClustersLay1 = 0;
96	fNClustersLay2 = 0;
97	fNTracklets = 0;
98
99	// loading the clusters
100	LoadClusterArrays(clusterTree);
101
102	// find the tracklets
103	AliDebug(1,"Looking for tracklets... ");
104
105	//###########################################################
106	// Loop on layer 1 : finding theta, phi and z
107	for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
108	Float_t x = fClustersLay1[iC1][0] - vtx[0];
109	Float_t y = fClustersLay1[iC1][1] - vtx[1];
110	Float_t z = fClustersLay1[iC1][2] - vtx[2];
111
112	Float_t r = TMath::Sqrt(TMath::Power(x,2) +
113	TMath::Power(y,2) +
114	TMath::Power(z,2));
115
116	fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
117	fClustersLay1[iC1][1] = TMath::ATan(y/x); // Store Phi
118	fClustersLay1[iC1][2] = z/r; // Store scaled z
119	}
120
121	// Loop on layer 2 : finding theta, phi and r
122	for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
123	Float_t x = fClustersLay2[iC2][0] - vtx[0];
124	Float_t y = fClustersLay2[iC2][1] - vtx[1];
125	Float_t z = fClustersLay2[iC2][2] - vtx[2];
126
127	Float_t r = TMath::Sqrt(TMath::Power(x,2) +
128	TMath::Power(y,2) +
129	TMath::Power(z,2));
130
131	fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
132	fClustersLay2[iC2][1] = TMath::ATan(y/x); // Store Phi
133	fClustersLay2[iC2][2] = z; // Store z
134
135	// this only needs to be initialized for the fNClustersLay2 first associations
136	fAssociationFlag[iC2] = kFALSE;
137	}
138
139	//###########################################################
140	// Loop on layer 1
141	for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
142
143	// reset of variables for multiple candidates
144	Int_t iC2WithBestPhi = 0; // reset
145	Float_t dPhimin = 100.; // just to put a huge number!
146
147	// Loop on layer 2
148	for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
149
150	// The following excludes double associations
151	if (!fAssociationFlag[iC2]) {
152
153	// find the difference in angles
154	Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
155	Float_t dPhi = fClustersLay2[iC2][1] - fClustersLay1[iC1][1];
156
157	// find the difference in z (between linear projection from layer 1
158	// and the actual point: Dzeta= z1/r1*r2 -z2)
159	Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
160	Float_t dZeta = fClustersLay2[iC1][2]*r2 - fClustersLay2[iC2][2];
161
162	if (fHistOn) {
163	fhClustersDPhi->Fill(dPhi);
164	fhClustersDTheta->Fill(dTheta);
165	fhClustersDZeta->Fill(dZeta);
166	fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
167	}
168	// make "elliptical" cut in Phi and Zeta!
169	Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
170	if (d>1) continue;
171
172	//look for the minimum distance in Phi: the minimum is in iC2WithBestPhi
173	if (TMath::Abs(dPhi) < dPhimin) {
174	dPhimin = TMath::Abs(dPhi);
175	iC2WithBestPhi = iC2;
176	}
177	}
178	} // end of loop over clusters in layer 2
179
180	if (dPhimin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
181
182	if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestPhi] = kTRUE; // flag the association
183
184	// store the tracklet
185
186	// use the average theta from the clusters in the two layers
187	fTracklets[fNTracklets][0] = 0.5*(fClustersLay1[iC1][0]+fClustersLay2[iC2WithBestPhi][0]);
188	// use the phi from the clusters in the first layer
189	fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
190	// Store the difference between phi1 and phi2
191	fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestPhi][1];
192	fNTracklets++;
193
194	AliDebug(1,Form(" Adding tracklet candidate %d (cluster %d of layer 1 and %d of layer 2)", fNTracklets, iC1));
195	}
196	} // end of loop over clusters in layer 1
197
198	AliDebug(1,Form("%d tracklets found", fNTracklets));
199	}
200
201	//____________________________________________________________________
202	void
203	AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
204	// This method
205	// - gets the clusters from the cluster tree
206	// - convert them into global coordinates
207	// - store them in the internal arrays
208
209	AliDebug(1,"Loading clusters ...");
210
211	fNClustersLay1 = 0;
212	fNClustersLay2 = 0;
213
214	TClonesArray* itsClusters = new TClonesArray("AliITSclusterV2");
215	TBranch* itsClusterBranch=itsClusterTree->GetBranch("Clusters");
216	itsClusterBranch->SetAddress(&itsClusters);
217
218	Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
219
220	// loop over the its subdetectors
221	for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
222
223	if (!itsClusterTree->GetEvent(iIts))
224	continue;
225
226	Int_t nClusters = itsClusters->GetEntriesFast();
227
228	// stuff needed to get the global coordinates
229	Double_t rot[9]; fGeometry->GetRotMatrix(iIts,rot);
230	Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det);
231	Float_t tx,ty,tz; fGeometry->GetTrans(lay,lad,det,tx,ty,tz);
232
233	// Below:
234	// "alpha" is the angle from the global X-axis to the
235	// local GEANT X'-axis ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
236	// "phi" is the angle from the global X-axis to the
237	// local cluster X"-axis
238
239	Double_t alpha = TMath::ATan2(rot[1],rot[0])+TMath::Pi();
240	Double_t itsPhi = TMath::Pi()/2+alpha;
241
242	if (lay==1) itsPhi+=TMath::Pi();
243	Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi);
244	Double_t r=txcp+tysp;
245
246	// loop over clusters
247	while(nClusters--) {
248	AliITSclusterV2* cluster = (AliITSclusterV2*)itsClusters->UncheckedAt(nClusters);
249
250	if (cluster->GetLayer()>1)
251	continue;
252
253	Float_t x = rcp - cluster->GetY()sp;
254	Float_t y = rsp + cluster->GetY()cp;
255	Float_t z = cluster->GetZ();
256
257	if (cluster->GetLayer()==0) {
258	fClustersLay1[fNClustersLay1][0] = x;
259	fClustersLay1[fNClustersLay1][1] = y;
260	fClustersLay1[fNClustersLay1][2] = z;
261	fNClustersLay1++;
262	}
263	if (cluster->GetLayer()==1) {
264	fClustersLay2[fNClustersLay2][0] = x;
265	fClustersLay2[fNClustersLay2][1] = y;
266	fClustersLay2[fNClustersLay2][2] = z;
267	fNClustersLay2++;
268	}
269
270	}// end of cluster loop
271	} // end of its "subdetector" loop
272
273	AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
274	}
275	//____________________________________________________________________
276	void
277	AliITSMultReconstructor::SaveHists() {
278
279	if (!fHistOn)
280	return;
281
282	cout << "Saving histograms" << endl;
283
284	fhClustersDPhi->Write();
285	fhClustersDTheta->Write();
286	fhClustersDZeta->Write();
287	fhDPhiVsDThetaAll->Write();
288	fhDPhiVsDThetaAcc->Write();
289	}