[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.
//
// 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 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.
// 
//  Author :  Tiziano Virgili 
//  
//  
//
//____________________________________________________________________

#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),
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).

  fGeometry =0;

  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.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),
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::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];

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