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