[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[5];
    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
      // if equal label in both clusters found this label is assigned
      // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
      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];
        fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3+label2];
      }
      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] = fClustersLay1[iC1][3];
        fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3];
      }

      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];
de4c520e	113	fClustersLay2[i] = new Float_t[6];
0939e22a	114	fTracklets[i] = new Float_t[5];
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
0939e22a	368	// if equal label in both clusters found this label is assigned
0939e22a	369	// if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
de4c520e	370	Int_t label1 = 0;
	371	Int_t label2 = 0;
	372	while (label2 < 3)
	373	{
	374	if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
	375	break;
de4c520e	376	label1++;
	377	if (label1 == 3)
	378	{
	379	label1 = 0;
	380	label2++;
	381	}
	382	}
	383
	384	if (label2 < 3)
	385	{
	386	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));
	387	fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
0939e22a	388	fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3+label2];
de4c520e	389	}
	390	else
	391	{
	392	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));
0939e22a	393	fTracklets[fNTracklets][3] = fClustersLay1[iC1][3];
0939e22a	394	fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3];
de4c520e	395	}
de4c520e	396
3ef75756	397	if (fHistOn) {
	398	Float_t eta=fTracklets[fNTracklets][0];
	399	eta= TMath::Tan(eta/2.);
	400	eta=-TMath::Log(eta);
	401	fhetaTracklets->Fill(eta);
	402	fhphiTracklets->Fill(fTracklets[fNTracklets][1]);
	403	}
ac903f1b	404
3ef75756	405	AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
	406	AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1,
	407	iC2WithBestDist));
	408	fNTracklets++;
ac903f1b	409	}
3ef75756	410
	411	// Delete the following else if you do not want to save Clusters!
	412
de4c520e	413	else { // This means that the cluster has not been associated
3ef75756	414
	415	// store the cluster
	416
968e8539	417	fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
968e8539	418	fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
de4c520e	419	AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
968e8539	420	fNSingleCluster, iC1));
968e8539	421	fNSingleCluster++;
3ef75756	422	}
3ef75756	423
ac903f1b	424	} // end of loop over clusters in layer 1
	425
	426	AliDebug(1,Form("%d tracklets found", fNTracklets));
	427	}
	428
	429	//____________________________________________________________________
	430	void
	431	AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
	432	// This method
	433	// - gets the clusters from the cluster tree
	434	// - convert them into global coordinates
	435	// - store them in the internal arrays
	436
	437	AliDebug(1,"Loading clusters ...");
	438
	439	fNClustersLay1 = 0;
	440	fNClustersLay2 = 0;
	441
b51872de	442	TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
b51872de	443	TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
ddced3c8	444
ac903f1b	445	itsClusterBranch->SetAddress(&itsClusters);
ddced3c8	446
ac903f1b	447	Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
f606f16a	448	Float_t cluGlo[3]={0.,0.,0.};
ddced3c8	449
ac903f1b	450	// loop over the its subdetectors
	451	for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
	452
	453	if (!itsClusterTree->GetEvent(iIts))
	454	continue;
	455
	456	Int_t nClusters = itsClusters->GetEntriesFast();
	457
ac903f1b	458	// loop over clusters
ac903f1b	459	while(nClusters--) {
de4c520e	460	AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
ac903f1b	461
	462	if (cluster->GetLayer()>1)
	463	continue;
	464
f606f16a	465	cluster->GetGlobalXYZ(cluGlo);
	466	Float_t x = cluGlo[0];
	467	Float_t y = cluGlo[1];
	468	Float_t z = cluGlo[2];
ac903f1b	469
	470	if (cluster->GetLayer()==0) {
	471	fClustersLay1[fNClustersLay1][0] = x;
	472	fClustersLay1[fNClustersLay1][1] = y;
	473	fClustersLay1[fNClustersLay1][2] = z;
de4c520e	474	for (Int_t i=0; i<3; i++)
de4c520e	475	fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
ac903f1b	476	fNClustersLay1++;
ac903f1b	477	}
de4c520e	478	if (cluster->GetLayer()==1) {
ac903f1b	479	fClustersLay2[fNClustersLay2][0] = x;
	480	fClustersLay2[fNClustersLay2][1] = y;
	481	fClustersLay2[fNClustersLay2][2] = z;
de4c520e	482	for (Int_t i=0; i<3; i++)
de4c520e	483	fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
ac903f1b	484	fNClustersLay2++;
	485	}
	486
	487	}// end of cluster loop
	488	} // end of its "subdetector" loop
	489
cedf398d	490	if (itsClusters) {
	491	itsClusters->Delete();
	492	delete itsClusters;
	493	itsClusters = 0;
	494	}
ac903f1b	495	AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2));
	496	}
	497	//____________________________________________________________________
	498	void
	499	AliITSMultReconstructor::SaveHists() {
3ef75756	500	// This method save the histograms on the output file
3ef75756	501	// (only if fHistOn is TRUE).
ac903f1b	502
	503	if (!fHistOn)
	504	return;
	505
ddced3c8	506	fhClustersDPhiAll->Write();
	507	fhClustersDThetaAll->Write();
	508	fhClustersDZetaAll->Write();
ac903f1b	509	fhDPhiVsDThetaAll->Write();
ddced3c8	510	fhDPhiVsDZetaAll->Write();
	511
	512	fhClustersDPhiAcc->Write();
	513	fhClustersDThetaAcc->Write();
	514	fhClustersDZetaAcc->Write();
ac903f1b	515	fhDPhiVsDThetaAcc->Write();
ddced3c8	516	fhDPhiVsDZetaAcc->Write();
	517
	518	fhetaTracklets->Write();
	519	fhphiTracklets->Write();
	520	fhetaClustersLay1->Write();
	521	fhphiClustersLay1->Write();
ac903f1b	522	}
de4c520e	523
eefb3acc	524