From ac903f1b3a3d49cedaca16dfc5e899181f4cca5e Mon Sep 17 00:00:00 2001 From: hristov Date: Wed, 5 Apr 2006 14:36:47 +0000 Subject: [PATCH] Reconstruction of particle multiplicity (T.Virgili, C.Jorgensen) --- ITS/AliITSMultReconstructor.cxx | 289 ++++++++++++++++++++++++++++++++ ITS/AliITSMultReconstructor.h | 87 ++++++++++ ITS/ITSrecLinkDef.h | 3 +- ITS/libITSrec.pkg | 3 +- ITS/testITSMultReco.C | 149 ++++++++++++++++ 5 files changed, 529 insertions(+), 2 deletions(-) create mode 100644 ITS/AliITSMultReconstructor.cxx create mode 100644 ITS/AliITSMultReconstructor.h create mode 100755 ITS/testITSMultReco.C diff --git a/ITS/AliITSMultReconstructor.cxx b/ITS/AliITSMultReconstructor.cxx new file mode 100644 index 00000000000..99cb81759a7 --- /dev/null +++ b/ITS/AliITSMultReconstructor.cxx @@ -0,0 +1,289 @@ +//____________________________________________________________________ +// +// 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. +// +// ----------------------------------------------------------------- +// +// TODO: +// +// - Introduce a rough pt estimation from the difference in phi ? +// - Allow for a more refined selection criterium in case of multiple +// candidates (for instance by introducing weights for the difference +// in Phi and Zeta). +// +//____________________________________________________________________ + +#include "AliITSMultReconstructor.h" + +#include "TTree.h" +#include "TH1F.h" +#include "TH2F.h" + + +#include "AliITSclusterV2.h" +#include "AliITSgeom.h" +#include "AliLog.h" + +//____________________________________________________________________ +ClassImp(AliITSMultReconstructor); + +//____________________________________________________________________ +AliITSMultReconstructor::AliITSMultReconstructor() { + + fGeometry =0; + + SetHistOn(); + SetPhiWindow(); + SetZetaWindow(); + SetOnlyOneTrackletPerC2(); + + fClustersLay1 = new Float_t*[300000]; + fClustersLay2 = new Float_t*[300000]; + fTracklets = new Float_t*[300000]; + fAssociationFlag = new Bool_t[300000]; + + for(Int_t i=0; i<300000; i++) { + fClustersLay1[i] = new Float_t[3]; + fClustersLay2[i] = new Float_t[3]; + fTracklets[i] = new Float_t[3]; + fAssociationFlag[i] = kFALSE; + } + + // definition of histograms + fhClustersDPhi = new TH1F("dphi", "dphi", 200,-0.1,0.1); + fhClustersDPhi->SetDirectory(0); + fhClustersDTheta = new TH1F("dtheta","dtheta",200,-0.1,0.1); + fhClustersDTheta->SetDirectory(0); + fhClustersDZeta = new TH1F("dzeta","dzeta",200,-0.2,0.2); + fhClustersDZeta->SetDirectory(0); + + fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",200,-0.1,0.1,200,-0.1,0.1); + fhDPhiVsDThetaAll->SetDirectory(0); + fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",200,-0.1,0.1,200,-0.1,0.1); + fhDPhiVsDThetaAcc->SetDirectory(0); + +} + + +//____________________________________________________________________ +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. + + cout << " HEEEEEEEEEEEEEEEEEEEEE " << flush << endl; + + // reset counters + fNClustersLay1 = 0; + fNClustersLay2 = 0; + fNTracklets = 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; iC1Fill(dPhi); + fhClustersDTheta->Fill(dTheta); + fhClustersDZeta->Fill(dZeta); + fhDPhiVsDThetaAll->Fill(dTheta, 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 in Phi: the minimum is in iC2WithBestPhi + if (TMath::Abs(dPhi) < dPhimin) { + dPhimin = TMath::Abs(dPhi); + iC2WithBestPhi = iC2; + } + } + } // end of loop over clusters in layer 2 + + if (dPhimin<100) { // This means that a cluster in layer 2 was found that mathes with iC1 + + if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestPhi] = kTRUE; // flag the association + + // store the tracklet + + // use the average theta from the clusters in the two layers + fTracklets[fNTracklets][0] = 0.5*(fClustersLay1[iC1][0]+fClustersLay2[iC2WithBestPhi][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[iC2WithBestPhi][1]; + fNTracklets++; + + AliDebug(1,Form(" Adding tracklet candidate %d (cluster %d of layer 1 and %d of layer 2)", fNTracklets, iC1)); + } + } // 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("AliITSclusterV2"); + TBranch* itsClusterBranch=itsClusterTree->GetBranch("Clusters"); + 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--) { + AliITSclusterV2* cluster = (AliITSclusterV2*)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; + fNClustersLay1++; + } + if (cluster->GetLayer()==1) { + fClustersLay2[fNClustersLay2][0] = x; + fClustersLay2[fNClustersLay2][1] = y; + fClustersLay2[fNClustersLay2][2] = z; + fNClustersLay2++; + } + + }// end of cluster loop + } // end of its "subdetector" loop + + AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2)); +} +//____________________________________________________________________ +void +AliITSMultReconstructor::SaveHists() { + + if (!fHistOn) + return; + + cout << "Saving histograms" << endl; + + fhClustersDPhi->Write(); + fhClustersDTheta->Write(); + fhClustersDZeta->Write(); + fhDPhiVsDThetaAll->Write(); + fhDPhiVsDThetaAcc->Write(); +} diff --git a/ITS/AliITSMultReconstructor.h b/ITS/AliITSMultReconstructor.h new file mode 100644 index 00000000000..ea2b1f313e2 --- /dev/null +++ b/ITS/AliITSMultReconstructor.h @@ -0,0 +1,87 @@ +#ifndef ALIITSMULTRECONSTRUCTOR_H +#define ALIITSMULTRECONSTRUCTOR_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $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 boolean fOnlyOneTrackletPerC2 allows to control if two clusters +// in layer 2 can be associated to the same cluster in layer 1 or not. +// +///////////////////////////////////////////////////////////////////////// + +#include "TObject.h" + +class TTree; +class TH1F; +class TH2F; +class AliITSgeom; + +class AliITSMultReconstructor : public TObject +{ +public: + AliITSMultReconstructor(); + + void SetGeometry(AliITSgeom* geo) {fGeometry = geo;} + + void Reconstruct(TTree* tree, Float_t* vtx, Float_t* vtxRes); + + void SetPhiWindow(Float_t w=0.08) {fPhiWindow=w;} + void SetZetaWindow(Float_t w=0.1) {fZetaWindow=w;} + void SetOnlyOneTrackletPerC2(Float_t b = kFALSE) {fOnlyOneTrackletPerC2 = b;} + + Int_t GetNClustersLayer1() const {return fNClustersLay1;} + Int_t GetNClustersLayer2() const {return fNClustersLay2;} + Int_t GetNTracklets() const {return fNTracklets;} + + Float_t* GetClusterLayer1(Int_t n) {return fClustersLay1[n];} + Float_t* GetClusterLayer2(Int_t n) {return fClustersLay2[n];} + Float_t* GetTracklet(Int_t n) {return fTracklets[n];} + + void SetHistOn(Bool_t b=kFALSE) {fHistOn=b;} + void SaveHists(); +protected: + + AliITSgeom* fGeometry; // ITS geometry + + Float_t** fClustersLay1; // clusters in the 1st layer of ITS + Float_t** fClustersLay2; // clusters in the 2nd layer of ITS + Float_t** fTracklets; // tracklets + Bool_t* fAssociationFlag; // flag for the associations + + Int_t fNClustersLay1; // Number of clusters (Layer1) + Int_t fNClustersLay2; // Number of clusters (Layer2) + Int_t fNTracklets; // Number of tracklets + + Float_t fPhiWindow; // Search window in phi + Float_t fZetaWindow; // SEarch window in eta + + Bool_t fOnlyOneTrackletPerC2; // only one tracklet per cluster in layer 2? + + Bool_t fHistOn; + + TH1F* fhClustersDPhi; + TH1F* fhClustersDTheta; + TH1F* fhClustersDZeta; + + TH2F* fhDPhiVsDThetaAll; + TH2F* fhDPhiVsDThetaAcc; + + void LoadClusterArrays(TTree* tree); + + ClassDef(AliITSMultReconstructor,0) +}; + +#endif diff --git a/ITS/ITSrecLinkDef.h b/ITS/ITSrecLinkDef.h index 2b9016713a7..dc5bafcc008 100644 --- a/ITS/ITSrecLinkDef.h +++ b/ITS/ITSrecLinkDef.h @@ -81,6 +81,7 @@ #pragma link C++ class AliITSBeamTestDigSDD+; #pragma link C++ class AliITSBeamTestDigSSD+; #pragma link C++ class AliITSBeamTestDigitizer+; - +//multiplicity +#pragma link C++ class AliITSMultReconstructor+; #endif diff --git a/ITS/libITSrec.pkg b/ITS/libITSrec.pkg index fad8f25ec3c..83bc77f3d36 100644 --- a/ITS/libITSrec.pkg +++ b/ITS/libITSrec.pkg @@ -49,7 +49,8 @@ SRCS = AliITSDetTypeRec.cxx \ AliITSBeamTestDigSDD.cxx \ AliITSBeamTestDigSPD.cxx \ AliITSBeamTestDigSSD.cxx \ - AliITSBeamTestDigitizer.cxx\ + AliITSBeamTestDigitizer.cxx \ + AliITSMultReconstructor.cxx HDRS:= $(SRCS:.cxx=.h) diff --git a/ITS/testITSMultReco.C b/ITS/testITSMultReco.C new file mode 100755 index 00000000000..49572f63c79 --- /dev/null +++ b/ITS/testITSMultReco.C @@ -0,0 +1,149 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) + +#include +#include +#include +#include + +#include "AliRunLoader.h" +#include "AliESD.h" +#include "AliRun.h" + +#include "AliITS.h" +#include "AliITSgeom.h" +#include "AliITSLoader.h" +#include "AliITSMultReconstructor.h" + +#endif + +void testITSMultReco(Char_t* dir = ".") { + + Char_t str[256]; + + // ######################################################## + // defining pointers + AliRunLoader* runLoader; + TFile* esdFile = 0; + TTree* esdTree = 0; + AliESD* esd = 0; + + // ######################################################### + // setup galice and runloader + + if (gAlice) { + delete gAlice->GetRunLoader(); + delete gAlice; + gAlice=0; + } + + sprintf(str,"%s/galice.root",dir); + runLoader = AliRunLoader::Open(str); + if (runLoader == 0x0) { + cout << "Can not open session"<LoadgAlice(); + + gAlice = runLoader->GetAliRun(); + runLoader->LoadKinematics(); + runLoader->LoadHeader(); + + // ######################################################### + // open esd file and get the tree + + // close it first to avoid memory leak + if (esdFile) + if (esdFile->IsOpen()) + esdFile->Close(); + + sprintf(str,"%s/AliESDs.root",dir); + esdFile = TFile::Open(str); + esdTree = (TTree*)esdFile->Get("esdTree"); + TBranch * esdBranch = esdTree->GetBranch("ESD"); + esdBranch->SetAddress(&esd); + + + // ######################################################### + // setup its stuff + + AliITS* its=(AliITS*)runLoader->GetAliRun()->GetDetector("ITS"); + if (!its) { + cout << " Can't get the ITS!" << endl; + return ; + } + AliITSgeom* itsGeo=its->GetITSgeom(); + if (!itsGeo) { + cout << " Can't get the ITS geometry!" << endl; + return ; + } + AliITSLoader* itsLoader = (AliITSLoader*)runLoader->GetLoader("ITSLoader"); + if (!itsLoader) { + cout << " Can't get the ITS loader!" << endl; + return ; + } + itsLoader->LoadRecPoints("read"); + + // ######################################################### + AliITSMultReconstructor* multReco = new AliITSMultReconstructor(); + multReco->SetGeometry(itsGeo); + + // ######################################################### + // getting number of events + + Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents(); + Int_t nESDEvents = esdBranch->GetEntries(); + + if (nEvents!=nESDEvents) { + cout << " Different number of events from runloader and esdtree!!!" + << nEvents << " / " << nESDEvents << endl; + return; + } + + // ######################################################## + // loop over number of events + cout << nEvents << " event(s) found in the file set" << endl; + for(Int_t i=0; iGetEvent(i); + esdBranch->GetEntry(i); + + // ######################################################## + // get the EDS vertex + const AliESDVertex* vtxESD = esd->GetVertex(); + Double_t vtx[3]; + vtxESD->GetXYZ(vtx); + Float_t esdVtx[3]; + esdVtx[0] = vtx[0]; + esdVtx[1] = vtx[1]; + esdVtx[2] = vtx[2]; + + ///######################################################### + // get ITS clusters + TTree* itsClusterTree = itsLoader->TreeR(); + if (!itsClusterTree) { + cerr<< " Can't get the ITS cluster tree !\n"; + return; + } + multReco->SetHistOn(kTRUE); + multReco->Reconstruct(itsClusterTree, esdVtx, esdVtx); + + + for (Int_t t=0; tGetNTracklets(); t++) { + + cout << " tracklet " << t + << " , theta = " << multReco->GetTracklet(t)[0] + << " , phi = " << multReco->GetTracklet(t)[1] << endl; + } + + } + + TFile* fout = new TFile("out.root","RECREATE"); + + multReco->SaveHists(); + fout->Write(); + fout->Close(); + + +} -- 2.39.3