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

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//____________________________________________________________________
// 
// AliITSMultReconstructor - find clusters in the pixels (theta and
// phi) and tracklets.
// 
// These can be used to extract charged particles multiplcicity from the ITS.
//
// A tracklet consist of two ITS clusters, one in the first pixel
// layer and one in the second. The clusters are associates if the 
// differencies in Phi (azimuth) and Zeta (longitudinal) are inside 
// a fiducial volume. In case of multiple candidates it is selected the
// candidate with minimum distance in Phi. 
// The parameter AssociationChoice allows to control if two clusters 
// in layer 2 can be associated to the same cluster in layer 1 or not.
//
// -----------------------------------------------------------------
// 
// 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 <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TTree.h>

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

//____________________________________________________________________
ClassImp(AliITSMultReconstructor)


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