[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 
//  
//  
//
//____________________________________________________________________

#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.,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,-3.14159,3.14159);
  fhphiTracklets->SetDirectory(0);
  fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
  fhetaClustersLay1->SetDirectory(0);
  fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100,-3.141,3.141);
  fhphiClustersLay1->SetDirectory(0);
}

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