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

/**************************************************************************
 * Copyright(c) 2007-2009, 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.
// (TRUE means double associations exluded; default = TRUE)
//
// Two methods return the number of traklets and the number of clusters 
// in the first SPD layer (GetNTracklets GetNSingleClusters)
//
// -----------------------------------------------------------------
// 
// NOTE: The cuts on phi and zeta depend 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.
// 
//  Author :  Tiziano Virgili 
//  
//  Recent updates (D. Elia, INFN Bari):
//     - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE) 
//     - phi definition changed to ALICE convention (0,2*TMath::pi()) 
//     - cluster coordinates taken with GetGlobalXYZ()
//
//  fGeometry removed
//____________________________________________________________________

#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():
fClustersLay1(0),
fClustersLay2(0),
fTracklets(0),
fSClusters(0),
fAssociationFlag(0),
fNClustersLay1(0),
fNClustersLay2(0),
fNTracklets(0),
fNSingleCluster(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).


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

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

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

  // definition of histograms
  fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,0.,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.,0.1);
  fhDPhiVsDZetaAcc->SetDirectory(0);
  fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
  fhDPhiVsDThetaAcc->SetDirectory(0);

  fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,0.0,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.,0.5);
  fhDPhiVsDZetaAll->SetDirectory(0);
  fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
  fhDPhiVsDThetaAll->SetDirectory(0);

  fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
  fhetaTracklets->SetDirectory(0);
  fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100, 0., 2*TMath::Pi());
  fhphiTracklets->SetDirectory(0);
  fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
  fhetaClustersLay1->SetDirectory(0);
  fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
  fhphiClustersLay1->SetDirectory(0);
}

//______________________________________________________________________
AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
fClustersLay1(mr.fClustersLay1),
fClustersLay2(mr.fClustersLay2),
fTracklets(mr.fTracklets),
fSClusters(mr.fSClusters),
fAssociationFlag(mr.fAssociationFlag),
fNClustersLay1(mr.fNClustersLay1),
fNClustersLay2(mr.fNClustersLay2),
fNTracklets(mr.fNTracklets),
fNSingleCluster(mr.fNSingleCluster),
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 [] fSClusters[i];
  }
  delete [] fClustersLay1;
  delete [] fClustersLay2;
  delete [] fTracklets;
  delete [] fSClusters;

  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; 
  fNSingleCluster = 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::Pi() + TMath::ATan2(-y,-x);  // 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::Pi() + TMath::ATan2(-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  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   = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
        // take into account boundary condition
        if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;	

	// 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];

      // find label
      Int_t label1 = 0;
      Int_t label2 = 0;
      while (label2 < 3)
      {
        if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
          break;
        label1++;
        if (label1 == 3)
        {
          label1 = 0;
          label2++;
        }
      }

      if (label2 < 3)
      {
        AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) fClustersLay1[iC1][3+label1], (Int_t) fClustersLay2[iC2WithBestDist][3+label2], fNTracklets));
        fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
      }
      else
      {
        AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", (Int_t) fClustersLay1[iC1][3], (Int_t) fClustersLay1[iC1][4], (Int_t) fClustersLay1[iC1][5], (Int_t) fClustersLay2[iC2WithBestDist][3], (Int_t) fClustersLay2[iC2WithBestDist][4], (Int_t) fClustersLay2[iC2WithBestDist][5], fNTracklets));
        fTracklets[fNTracklets][3] = -2;
      }

      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
       
      fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
      fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
      AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
		      fNSingleCluster, iC1));
      fNSingleCluster++;
    }

  } // 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();  
  Float_t cluGlo[3]={0.,0.,0.};
 
  // loop over the its subdetectors
  for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
    
    if (!itsClusterTree->GetEvent(iIts)) 
      continue;
    
    Int_t nClusters = itsClusters->GetEntriesFast();
    
    // loop over clusters
    while(nClusters--) {
      AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
      
      if (cluster->GetLayer()>1) 
	continue;            
      
      cluster->GetGlobalXYZ(cluGlo);
      Float_t x = cluGlo[0];
      Float_t y = cluGlo[1];
      Float_t z = cluGlo[2];      
      
      if (cluster->GetLayer()==0) {
	fClustersLay1[fNClustersLay1][0] = x;
	fClustersLay1[fNClustersLay1][1] = y;
	fClustersLay1[fNClustersLay1][2] = z;
	for (Int_t i=0; i<3; i++)
		fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
	fNClustersLay1++;
      }
      if (cluster->GetLayer()==1) {
	fClustersLay2[fNClustersLay2][0] = x;
	fClustersLay2[fNClustersLay2][1] = y;
	fClustersLay2[fNClustersLay2][2] = z;
	for (Int_t i=0; i<3; i++)
		fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
	fNClustersLay2++;
      }
      
    }// end of cluster loop
  } // end of its "subdetector" loop  
  
  if (itsClusters) {
    itsClusters->Delete();
    delete itsClusters;
    itsClusters = 0;
  }
  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	/**************************************************************************
eefb3acc	2	* Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
7ca4655f	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.
de4c520e	30	// The parameter AssociationChoice allows to control if two clusters
ac903f1b	31	// in layer 2 can be associated to the same cluster in layer 1 or not.
02a95988	32	// (TRUE means double associations exluded; default = TRUE)
ac903f1b	33	//
968e8539	34	// Two methods return the number of traklets and the number of clusters
	35	// in the first SPD layer (GetNTracklets GetNSingleClusters)
	36	//
ac903f1b	37	// -----------------------------------------------------------------
ac903f1b	38	//
02a95988	39	// NOTE: The cuts on phi and zeta depend on the interacting system (p-p
3ef75756	40	// or Pb-Pb). Please, check the file AliITSMultReconstructor.h and be
3ef75756	41	// sure that SetPhiWindow and SetZetaWindow are defined accordingly.
ac903f1b	42	//
968e8539	43	// Author : Tiziano Virgili
3ef75756	44	//
f606f16a	45	// Recent updates (D. Elia, INFN Bari):
	46	// - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE)
	47	// - phi definition changed to ALICE convention (0,2*TMath::pi())
	48	// - cluster coordinates taken with GetGlobalXYZ()
ac903f1b	49	//
f606f16a	50	// fGeometry removed
ac903f1b	51	//____________________________________________________________________
ac903f1b	52
7ca4655f	53	#include <TClonesArray.h>
	54	#include <TH1F.h>
	55	#include <TH2F.h>
	56	#include <TTree.h>
ac903f1b	57
7ca4655f	58	#include "AliITSMultReconstructor.h"
b51872de	59	#include "AliITSRecPoint.h"
ac903f1b	60	#include "AliITSgeom.h"
	61	#include "AliLog.h"
	62
	63	//____________________________________________________________________
0762f3a8	64	ClassImp(AliITSMultReconstructor)
ac903f1b	65
3ef75756	66
ac903f1b	67	//____________________________________________________________________
7537d03c	68	AliITSMultReconstructor::AliITSMultReconstructor():
7537d03c	69	fClustersLay1(0),
	70	fClustersLay2(0),
	71	fTracklets(0),
968e8539	72	fSClusters(0),
7537d03c	73	fAssociationFlag(0),
	74	fNClustersLay1(0),
	75	fNClustersLay2(0),
	76	fNTracklets(0),
968e8539	77	fNSingleCluster(0),
7537d03c	78	fPhiWindow(0),
	79	fZetaWindow(0),
	80	fOnlyOneTrackletPerC2(0),
	81	fHistOn(0),
	82	fhClustersDPhiAcc(0),
	83	fhClustersDThetaAcc(0),
	84	fhClustersDZetaAcc(0),
	85	fhClustersDPhiAll(0),
	86	fhClustersDThetaAll(0),
	87	fhClustersDZetaAll(0),
	88	fhDPhiVsDThetaAll(0),
	89	fhDPhiVsDThetaAcc(0),
	90	fhDPhiVsDZetaAll(0),
	91	fhDPhiVsDZetaAcc(0),
	92	fhetaTracklets(0),
	93	fhphiTracklets(0),
	94	fhetaClustersLay1(0),
	95	fhphiClustersLay1(0){
3ef75756	96	// Method to reconstruct the charged particles multiplicity with the
3ef75756	97	// SPD (tracklets).
ac903f1b	98
ac903f1b	99
	100	SetHistOn();
	101	SetPhiWindow();
	102	SetZetaWindow();
	103	SetOnlyOneTrackletPerC2();
	104
	105	fClustersLay1 = new Float_t*[300000];
	106	fClustersLay2 = new Float_t*[300000];
	107	fTracklets = new Float_t*[300000];
968e8539	108	fSClusters = new Float_t*[300000];
ac903f1b	109	fAssociationFlag = new Bool_t[300000];
	110
	111	for(Int_t i=0; i<300000; i++) {
de4c520e	112	fClustersLay1[i] = new Float_t[6];
	113	fClustersLay2[i] = new Float_t[6];
	114	fTracklets[i] = new Float_t[4];
968e8539	115	fSClusters[i] = new Float_t[2];
ac903f1b	116	fAssociationFlag[i] = kFALSE;
	117	}
	118
	119	// definition of histograms
02a95988	120	fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1);
ddced3c8	121	fhClustersDPhiAcc->SetDirectory(0);
	122	fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
	123	fhClustersDThetaAcc->SetDirectory(0);
	124	fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
	125	fhClustersDZetaAcc->SetDirectory(0);
	126
02a95988	127	fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
ddced3c8	128	fhDPhiVsDZetaAcc->SetDirectory(0);
02a95988	129	fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
ac903f1b	130	fhDPhiVsDThetaAcc->SetDirectory(0);
ac903f1b	131
02a95988	132	fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5);
ddced3c8	133	fhClustersDPhiAll->SetDirectory(0);
	134	fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
	135	fhClustersDThetaAll->SetDirectory(0);
	136	fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
	137	fhClustersDZetaAll->SetDirectory(0);
	138
02a95988	139	fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
ddced3c8	140	fhDPhiVsDZetaAll->SetDirectory(0);
02a95988	141	fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
ddced3c8	142	fhDPhiVsDThetaAll->SetDirectory(0);
	143
	144	fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.);
96c2c35d	145	fhetaTracklets->SetDirectory(0);
f606f16a	146	fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi());
96c2c35d	147	fhphiTracklets->SetDirectory(0);
ddced3c8	148	fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.);
96c2c35d	149	fhetaClustersLay1->SetDirectory(0);
f606f16a	150	fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
96c2c35d	151	fhphiClustersLay1->SetDirectory(0);
ac903f1b	152	}
ddced3c8	153
3ef75756	154	//______________________________________________________________________
7537d03c	155	AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
7537d03c	156	fClustersLay1(mr.fClustersLay1),
	157	fClustersLay2(mr.fClustersLay2),
	158	fTracklets(mr.fTracklets),
968e8539	159	fSClusters(mr.fSClusters),
7537d03c	160	fAssociationFlag(mr.fAssociationFlag),
	161	fNClustersLay1(mr.fNClustersLay1),
	162	fNClustersLay2(mr.fNClustersLay2),
	163	fNTracklets(mr.fNTracklets),
968e8539	164	fNSingleCluster(mr.fNSingleCluster),
7537d03c	165	fPhiWindow(mr.fPhiWindow),
	166	fZetaWindow(mr.fZetaWindow),
	167	fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
	168	fHistOn(mr.fHistOn),
	169	fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
	170	fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
	171	fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
	172	fhClustersDPhiAll(mr.fhClustersDPhiAll),
	173	fhClustersDThetaAll(mr.fhClustersDThetaAll),
	174	fhClustersDZetaAll(mr.fhClustersDZetaAll),
	175	fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
	176	fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
	177	fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
	178	fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
	179	fhetaTracklets(mr.fhetaTracklets),
	180	fhphiTracklets(mr.fhphiTracklets),
	181	fhetaClustersLay1(mr.fhetaClustersLay1),
	182	fhphiClustersLay1(mr.fhphiClustersLay1) {
3ef75756	183	// Copy constructor
7537d03c	184
3ef75756	185	}
	186
	187	//______________________________________________________________________
7537d03c	188	AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
3ef75756	189	// Assignment operator
7537d03c	190	this->~AliITSMultReconstructor();
7537d03c	191	new(this) AliITSMultReconstructor(mr);
3ef75756	192	return *this;
	193	}
	194
	195	//______________________________________________________________________
	196	AliITSMultReconstructor::~AliITSMultReconstructor(){
	197	// Destructor
1ba5b31c	198
	199	// delete histograms
	200	delete fhClustersDPhiAcc;
	201	delete fhClustersDThetaAcc;
	202	delete fhClustersDZetaAcc;
	203	delete fhClustersDPhiAll;
	204	delete fhClustersDThetaAll;
	205	delete fhClustersDZetaAll;
	206	delete fhDPhiVsDThetaAll;
	207	delete fhDPhiVsDThetaAcc;
	208	delete fhDPhiVsDZetaAll;
	209	delete fhDPhiVsDZetaAcc;
	210	delete fhetaTracklets;
	211	delete fhphiTracklets;
	212	delete fhetaClustersLay1;
	213	delete fhphiClustersLay1;
	214
	215	// delete arrays
	216	for(Int_t i=0; i<300000; i++) {
	217	delete [] fClustersLay1[i];
	218	delete [] fClustersLay2[i];
	219	delete [] fTracklets[i];
968e8539	220	delete [] fSClusters[i];
ddced3c8	221	}
1ba5b31c	222	delete [] fClustersLay1;
	223	delete [] fClustersLay2;
	224	delete [] fTracklets;
968e8539	225	delete [] fSClusters;
1ba5b31c	226
1ba5b31c	227	delete [] fAssociationFlag;
ddced3c8	228	}
ac903f1b	229
	230	//____________________________________________________________________
	231	void
	232	AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
	233	//
	234	// - calls LoadClusterArray that finds the position of the clusters
	235	// (in global coord)
	236	// - convert the cluster coordinates to theta, phi (seen from the
	237	// interaction vertex). The third coordinate is used for ....
	238	// - makes an array of tracklets
	239	//
	240	// After this method has been called, the clusters of the two layers
	241	// and the tracklets can be retrieved by calling the Get'er methods.
	242
ac903f1b	243	// reset counters
	244	fNClustersLay1 = 0;
	245	fNClustersLay2 = 0;
	246	fNTracklets = 0;
968e8539	247	fNSingleCluster = 0;
ac903f1b	248	// loading the clusters
ac903f1b	249	LoadClusterArrays(clusterTree);
3ef75756	250
ac903f1b	251	// find the tracklets
	252	AliDebug(1,"Looking for tracklets... ");
	253
	254	//###########################################################
	255	// Loop on layer 1 : finding theta, phi and z
	256	for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
	257	Float_t x = fClustersLay1[iC1][0] - vtx[0];
	258	Float_t y = fClustersLay1[iC1][1] - vtx[1];
	259	Float_t z = fClustersLay1[iC1][2] - vtx[2];
ddced3c8	260
ac903f1b	261	Float_t r = TMath::Sqrt(TMath::Power(x,2) +
	262	TMath::Power(y,2) +
	263	TMath::Power(z,2));
	264
eefb3acc	265	fClustersLay1[iC1][0] = TMath::ACos(z/r); // Store Theta
	266	fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
	267	fClustersLay1[iC1][2] = z/r; // Store scaled z
ddced3c8	268	if (fHistOn) {
	269	Float_t eta=fClustersLay1[iC1][0];
	270	eta= TMath::Tan(eta/2.);
	271	eta=-TMath::Log(eta);
	272	fhetaClustersLay1->Fill(eta);
de4c520e	273	fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
ddced3c8	274	}
96c2c35d	275	}
ac903f1b	276
	277	// Loop on layer 2 : finding theta, phi and r
	278	for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
	279	Float_t x = fClustersLay2[iC2][0] - vtx[0];
	280	Float_t y = fClustersLay2[iC2][1] - vtx[1];
	281	Float_t z = fClustersLay2[iC2][2] - vtx[2];
ddced3c8	282
ac903f1b	283	Float_t r = TMath::Sqrt(TMath::Power(x,2) +
	284	TMath::Power(y,2) +
	285	TMath::Power(z,2));
	286
eefb3acc	287	fClustersLay2[iC2][0] = TMath::ACos(z/r); // Store Theta
	288	fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x); // Store Phi
	289	fClustersLay2[iC2][2] = z; // Store z
ac903f1b	290
ddced3c8	291	// this only needs to be initialized for the fNClustersLay2 first associations
ac903f1b	292	fAssociationFlag[iC2] = kFALSE;
	293	}
	294
	295	//###########################################################
	296	// Loop on layer 1
	297	for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
	298
	299	// reset of variables for multiple candidates
ddced3c8	300	Int_t iC2WithBestDist = 0; // reset
3ef75756	301	Float_t distmin = 100.; // just to put a huge number!
ddced3c8	302	Float_t dPhimin = 0.; // Used for histograms only!
	303	Float_t dThetamin = 0.; // Used for histograms only!
	304	Float_t dZetamin = 0.; // Used for histograms only!
ac903f1b	305
	306	// Loop on layer 2
	307	for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
	308
	309	// The following excludes double associations
	310	if (!fAssociationFlag[iC2]) {
	311
	312	// find the difference in angles
	313	Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
02a95988	314	Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
	315	// take into account boundary condition
	316	if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;
	317
ac903f1b	318	// find the difference in z (between linear projection from layer 1
ac903f1b	319	// and the actual point: Dzeta= z1/r1*r2 -z2)
ddced3c8	320	Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
de4c520e	321	Float_t dZeta = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
ddced3c8	322
	323	if (fHistOn) {
	324	fhClustersDPhiAll->Fill(dPhi);
	325	fhClustersDThetaAll->Fill(dTheta);
	326	fhClustersDZetaAll->Fill(dZeta);
ac903f1b	327	fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
ddced3c8	328	fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
ac903f1b	329	}
	330	// make "elliptical" cut in Phi and Zeta!
	331	Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
3ef75756	332
ac903f1b	333	if (d>1) continue;
ac903f1b	334
ddced3c8	335	//look for the minimum distance: the minimum is in iC2WithBestDist
3ef75756	336	if (TMath::Sqrt(dZetadZeta+(r2dPhir2dPhi)) < distmin ) {
3ef75756	337	distmin=TMath::Sqrt(dZetadZeta + (r2dPhir2dPhi));
ddced3c8	338	dPhimin = dPhi;
	339	dThetamin = dTheta;
	340	dZetamin = dZeta;
	341	iC2WithBestDist = iC2;
ac903f1b	342	}
	343	}
	344	} // end of loop over clusters in layer 2
	345
3ef75756	346	if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
	347
	348	if (fHistOn) {
de4c520e	349	fhClustersDPhiAcc->Fill(dPhimin);
3ef75756	350	fhClustersDThetaAcc->Fill(dThetamin);
	351	fhClustersDZetaAcc->Fill(dZetamin);
	352	fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
	353	fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
	354	}
ac903f1b	355
ddced3c8	356	if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
ac903f1b	357
	358	// store the tracklet
	359
de4c520e	360	// use the theta from the clusters in the first layer
ddced3c8	361	fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
de4c520e	362	// use the phi from the clusters in the first layer
ac903f1b	363	fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
ac903f1b	364	// Store the difference between phi1 and phi2
de4c520e	365	fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
	366
	367	// find label
	368	Int_t label1 = 0;
	369	Int_t label2 = 0;
	370	while (label2 < 3)
	371	{
	372	if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
	373	break;
de4c520e	374	label1++;
	375	if (label1 == 3)
	376	{
	377	label1 = 0;
	378	label2++;
	379	}
	380	}
	381
	382	if (label2 < 3)
	383	{
	384	AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) fClustersLay1[iC1][3+label1], (Int_t) fClustersLay2[iC2WithBestDist][3+label2], fNTracklets));
	385	fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
	386	}
	387	else
	388	{
	389	AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", (Int_t) fClustersLay1[iC1][3], (Int_t) fClustersLay1[iC1][4], (Int_t) fClustersLay1[iC1][5], (Int_t) fClustersLay2[iC2WithBestDist][3], (Int_t) fClustersLay2[iC2WithBestDist][4], (Int_t) fClustersLay2[iC2WithBestDist][5], fNTracklets));
	390	fTracklets[fNTracklets][3] = -2;
	391	}
	392
3ef75756	393	if (fHistOn) {
	394	Float_t eta=fTracklets[fNTracklets][0];
	395	eta= TMath::Tan(eta/2.);
	396	eta=-TMath::Log(eta);
	397	fhetaTracklets->Fill(eta);
	398	fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
	399	}
ac903f1b	400
3ef75756	401	AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
	402	AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1,
	403	iC2WithBestDist));
	404	fNTracklets++;
ac903f1b	405	}
3ef75756	406
	407	// Delete the following else if you do not want to save Clusters!
	408
de4c520e	409	else { // This means that the cluster has not been associated
3ef75756	410
	411	// store the cluster
	412
968e8539	413	fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
968e8539	414	fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
de4c520e	415	AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
968e8539	416	fNSingleCluster, iC1));
968e8539	417	fNSingleCluster++;
3ef75756	418	}
3ef75756	419
ac903f1b	420	} // end of loop over clusters in layer 1
	421
	422	AliDebug(1,Form("%d tracklets found", fNTracklets));
	423	}
	424
	425	//____________________________________________________________________
	426	void
	427	AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
	428	// This method
	429	// - gets the clusters from the cluster tree
	430	// - convert them into global coordinates
	431	// - store them in the internal arrays
	432
	433	AliDebug(1,"Loading clusters ...");
	434
	435	fNClustersLay1 = 0;
	436	fNClustersLay2 = 0;
	437
b51872de	438	TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
b51872de	439	TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
ddced3c8	440
ac903f1b	441	itsClusterBranch->SetAddress(&itsClusters);
ddced3c8	442
ac903f1b	443	Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
f606f16a	444	Float_t cluGlo[3]={0.,0.,0.};
ddced3c8	445
ac903f1b	446	// loop over the its subdetectors
	447	for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
	448
	449	if (!itsClusterTree->GetEvent(iIts))
	450	continue;
	451
	452	Int_t nClusters = itsClusters->GetEntriesFast();
	453
ac903f1b	454	// loop over clusters
ac903f1b	455	while(nClusters--) {
de4c520e	456	AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
ac903f1b	457
	458	if (cluster->GetLayer()>1)
	459	continue;
	460
f606f16a	461	cluster->GetGlobalXYZ(cluGlo);
	462	Float_t x = cluGlo[0];
	463	Float_t y = cluGlo[1];
	464	Float_t z = cluGlo[2];
ac903f1b	465
	466	if (cluster->GetLayer()==0) {
	467	fClustersLay1[fNClustersLay1][0] = x;
	468	fClustersLay1[fNClustersLay1][1] = y;
	469	fClustersLay1[fNClustersLay1][2] = z;
de4c520e	470	for (Int_t i=0; i<3; i++)
de4c520e	471	fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
ac903f1b	472	fNClustersLay1++;
ac903f1b	473	}
de4c520e	474	if (cluster->GetLayer()==1) {
ac903f1b	475	fClustersLay2[fNClustersLay2][0] = x;
	476	fClustersLay2[fNClustersLay2][1] = y;
	477	fClustersLay2[fNClustersLay2][2] = z;
de4c520e	478	for (Int_t i=0; i<3; i++)
de4c520e	479	fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
ac903f1b	480	fNClustersLay2++;
	481	}
	482
	483	}// end of cluster loop
	484	} // end of its "subdetector" loop
	485
cedf398d	486	if (itsClusters) {
	487	itsClusters->Delete();
	488	delete itsClusters;
	489	itsClusters = 0;
	490	}
ac903f1b	491	AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
	492	}
	493	//____________________________________________________________________
	494	void
	495	AliITSMultReconstructor::SaveHists() {
3ef75756	496	// This method save the histograms on the output file
3ef75756	497	// (only if fHistOn is TRUE).
ac903f1b	498
	499	if (!fHistOn)
	500	return;
	501
ddced3c8	502	fhClustersDPhiAll->Write();
	503	fhClustersDThetaAll->Write();
	504	fhClustersDZetaAll->Write();
ac903f1b	505	fhDPhiVsDThetaAll->Write();
ddced3c8	506	fhDPhiVsDZetaAll->Write();
	507
	508	fhClustersDPhiAcc->Write();
	509	fhClustersDThetaAcc->Write();
	510	fhClustersDZetaAcc->Write();
ac903f1b	511	fhDPhiVsDThetaAcc->Write();
ddced3c8	512	fhDPhiVsDZetaAcc->Write();
	513
	514	fhetaTracklets->Write();
	515	fhphiTracklets->Write();
	516	fhetaClustersLay1->Write();
	517	fhphiClustersLay1->Write();
ac903f1b	518	}
de4c520e	519
eefb3acc	520