From 409f8c84935db0e639e3a3fb8a4447675ab1c326 Mon Sep 17 00:00:00 2001 From: barbera Date: Thu, 31 May 2001 21:02:56 +0000 Subject: [PATCH] Dubna model of SPD simulation. It was the default before --- ITS/AliITSClusterFinderSPDdubna.cxx | 402 ++++++++++++++++++ ITS/AliITSClusterFinderSPDdubna.h | 69 +++ ITS/AliITSFindClustersDubna.C | 101 +++++ ITS/AliITSHits2DigitsDubna.C | 109 +++++ ITS/AliITSresponseSPDdubna.cxx | 32 ++ ITS/AliITSresponseSPDdubna.h | 67 +++ ITS/AliITSsimulationSPDdubna.cxx | 625 ++++++++++++++++++++++++++++ ITS/AliITSsimulationSPDdubna.h | 51 +++ ITS/AliITStestDubna.C | 32 ++ ITS/ITSDigitsToClustersDubna.C | 161 +++++++ ITS/ITSHitsToDigitsDubna.C | 177 ++++++++ ITS/ITSreadClustTestSPDDubna.C | 111 +++++ ITS/SPDclusterTestDubna.C | 565 +++++++++++++++++++++++++ 13 files changed, 2502 insertions(+) create mode 100644 ITS/AliITSClusterFinderSPDdubna.cxx create mode 100644 ITS/AliITSClusterFinderSPDdubna.h create mode 100644 ITS/AliITSFindClustersDubna.C create mode 100644 ITS/AliITSHits2DigitsDubna.C create mode 100644 ITS/AliITSresponseSPDdubna.cxx create mode 100644 ITS/AliITSresponseSPDdubna.h create mode 100644 ITS/AliITSsimulationSPDdubna.cxx create mode 100644 ITS/AliITSsimulationSPDdubna.h create mode 100644 ITS/AliITStestDubna.C create mode 100644 ITS/ITSDigitsToClustersDubna.C create mode 100644 ITS/ITSHitsToDigitsDubna.C create mode 100644 ITS/ITSreadClustTestSPDDubna.C create mode 100644 ITS/SPDclusterTestDubna.C diff --git a/ITS/AliITSClusterFinderSPDdubna.cxx b/ITS/AliITSClusterFinderSPDdubna.cxx new file mode 100644 index 00000000000..9d3d99c3d55 --- /dev/null +++ b/ITS/AliITSClusterFinderSPDdubna.cxx @@ -0,0 +1,402 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* +$Log$ +Revision 1.15 2001/05/01 22:37:44 nilsen +Added a dummy argument to FindRawClusters. Argument used for SSD version. + +Revision 1.14 2001/03/05 14:48:46 nilsen +Fixed a reoccuring bug. Array sizes must be declare const. + +*/ + +#include +#include "AliITSClusterFinderSPDdubna.h" +#include "AliITSMapA1.h" +#include "AliITS.h" +#include "AliITSdigit.h" +#include "AliITSRawCluster.h" +#include "AliITSRecPoint.h" +#include "AliITSsegmentation.h" +#include "AliITSresponse.h" +#include "AliRun.h" + + + +ClassImp(AliITSClusterFinderSPDdubna) + +//---------------------------------------------------------- +AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna +(AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp) +{ + // constructor + fSegmentation=seg; + fDigits=digits; + fClusters=recp; + fNclusters= 0; + fMap=new AliITSMapA1(fSegmentation,fDigits); + SetDx(); + SetDz(); + SetNCells(); +} + +//_____________________________________________________________________________ +AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna() +{ + // constructor + fSegmentation=0; + fDigits=0; + fClusters=0; + fNclusters=0; + fMap=0; + SetDx(); + SetDz(); + SetNCells(); + +} + +//_____________________________________________________________________________ +AliITSClusterFinderSPDdubna::~AliITSClusterFinderSPDdubna() +{ + // destructor + if (fMap) delete fMap; + + +} +//_____________________________________________________________________________ + +void AliITSClusterFinderSPDdubna::Find1DClusters() +{ + // Find one dimensional clusters, i.e. + // in r*phi(x) direction for each colunm in z direction + + AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); + + // retrieve the parameters + Int_t fNofPixels = fSegmentation->Npz(); + Int_t fMaxNofSamples = fSegmentation->Npx(); + + // read in digits -> do not apply threshold + // signal in fired pixels is always 1 + + fMap->FillMap(); + + Int_t nofFoundClusters = 0; + + Int_t k,it,m; + for(k=0;k= fMaxNofSamples) break; // ! no possible for the fadc + + if(fMap->TestHit(k,id) == kUnused) { // start of the cluster + lclx += 1; + if(lclx == 1) xstart = id; + + } + + if(lclx > 0 && fMap->TestHit(k,id) == kEmpty) { + // end of cluster if a gap exists + xstop = id-1; + ilcl = 1; + break; + } + + } // end of m-loop + + if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window, + // continue the "it" loop + + if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished + + if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) { + // cluster end is outside of the window, + mmax += 5; // increase mmax and repeat the cluster + // finding + it -= 1; + } + + if(id >= fMaxNofSamples && lclx > 0) { // the x row finished but + xstop = fMaxNofSamples - 1; // the end cluster exists + ilcl = 1; + } + + // --- Calculate z and x coordinates for one dimensional clusters + + if(ilcl == 1) { // new cluster exists + it = id; + mmax = 10; + nofFoundClusters++; + Float_t clusterCharge = 0.; + Float_t zpitch = fSegmentation->Dpz(k+1); + Float_t clusterZ, dummyX; + Int_t dummy=0; + fSegmentation->GetPadCxz(dummy,k,dummyX,clusterZ); + Float_t zstart = clusterZ - 0.5*zpitch; + Float_t zstop = clusterZ + 0.5*zpitch; + Float_t clusterX = 0.; + Int_t xstartfull = xstart; + Int_t xstopfull = xstop; + Int_t clusterSizeX = lclx; + Int_t clusterSizeZ = 1; + + Int_t its; + for(its=xstart; its<=xstop; its++) { + Int_t firedpixel=0; + if (fMap->GetHitIndex(k,its)>=0) firedpixel=1; + clusterCharge += firedpixel; + clusterX +=its + 0.5; + } + Float_t fRphiPitch = fSegmentation->Dpx(dummy); + clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x + + // Write the points (coordinates and some cluster information) to the + // AliITSRawClusterSPD object + + AliITSRawClusterSPD clust(clusterZ,clusterX,clusterCharge,clusterSizeZ,clusterSizeX,xstart,xstop,xstartfull,xstopfull,zstart,zstop,k); + + iTS->AddCluster(0,&clust); + + } // new cluster (ilcl=1) + } // X direction loop (it) + } // Z direction loop (k) + + //fMap->ClearMap(); + return; + +} + +//_____________________________________________________________________________ +void AliITSClusterFinderSPDdubna::GroupClusters() +{ + // Find two dimensional clusters, i.e. group one dimensional clusters + // into two dimensional ones (go both in x and z directions). + + // get number of clusters for this module + Int_t nofClusters = fClusters->GetEntriesFast(); + nofClusters -= fNclusters; + + AliITSRawClusterSPD *clusterI; + AliITSRawClusterSPD *clusterJ; + + Int_t *label=new Int_t[nofClusters]; + Int_t i,j; + for(i=0; iAt(i); + clusterJ = (AliITSRawClusterSPD*) fClusters->At(j); + Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx); + if(pair) { + + // if((clusterI->XStop() == clusterJ->XStart()-1)||(clusterI->XStart()==clusterJ->XStop()+1)) cout<<"!! Diagonal cluster"<PrintInfo(); + clusterJ->PrintInfo(); + */ + + clusterI->Add(clusterJ); + // cout << "remove cluster " << j << endl; + label[j] = 1; + fClusters->RemoveAt(j); + + /* + cout << "cluster " << i << " after grouping" << endl; + clusterI->PrintInfo(); + */ + + } // pair + } // J clusters + label[i] = 1; + } // I clusters + fClusters->Compress(); + + /* + Int_t totalNofClusters = fClusters->GetEntriesFast(); + cout << " Nomber of clusters at the group end ="<< totalNofClusters<GetEntriesFast(); + nofClusters -= fNclusters; + + Int_t i, ix, iz, jx, jz, xstart, xstop, zstart, zstop, nclx, nclz; + const Int_t trmax = 100; + Int_t cltracks[trmax], itr, tracki, ii, is, js, ie, ntr, tr0, tr1, tr2; + + for(i=0; iAt(i); + + nclx = clusterI->NclX(); + nclz = clusterI->NclZ(); + xstart = clusterI->XStartf(); + xstop = clusterI->XStopf(); + zstart = clusterI->Zend()-nclz+1; + zstop = clusterI->Zend(); + Int_t ind; + + for(iz=0; izGetHitIndex(jz,jx); + if(ind < 0) { + continue; + } + if(ind == 0 && iz >= 0 && ix > 0) { + continue; + } + if(ind == 0 && iz > 0 && ix >= 0) { + continue; + } + if(ind == 0 && iz == 0 && ix == 0 && i > 0) { + continue; + } + + AliITSdigitSPD *dig = (AliITSdigitSPD*)fMap->GetHit(jz,jx); + + /* + signal=dig->fSignal; + track0=dig->fTracks[0]; + track1=dig->fTracks[1]; + track2=dig->fTracks[2]; + */ + + for(itr=0; itr<3; itr++) { + tracki = dig->fTracks[itr]; + if(tracki >= 0) { + ii += 1; + if(ii > 90) { + } + if(ii < 99) cltracks[ii-1] = tracki; + } + } + } // ix pixel + } // iz pixel + + for(is=0; is= 0) { + ntr=ntr+1; + if(ntr==1) tr0=cltracks[ie]; + if(ntr==2) tr1=cltracks[ie]; + if(ntr==3) tr2=cltracks[ie]; + } + } + // if delta ray only + if(ntr == 0) ntr = 1; + + clusterI->SetNTracks(ntr); + clusterI->SetTracks(tr0,tr1,tr2); + + } // I cluster + +} +//_____________________________________________________________________________ + +void AliITSClusterFinderSPDdubna::GetRecPoints() +{ + // get rec points + AliITS *iTS=(AliITS*)gAlice->GetModule("ITS"); + + // get number of clusters for this module + Int_t nofClusters = fClusters->GetEntriesFast(); + nofClusters -= fNclusters; + const Float_t kconv = 1.0e-4; + const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3 + const Float_t kRMSz = 70.0*kconv; // microns -> cm ITS TDR Table 1.3 + + Float_t spdLength = fSegmentation->Dz(); + Float_t spdWidth = fSegmentation->Dx(); + + Int_t i; + Int_t track0, track1, track2; + + for(i=0; iAt(i); + clusterI->GetTracks(track0, track1, track2); + AliITSRecPoint rnew; + + rnew.SetX((clusterI->X() - spdWidth/2)*kconv); + rnew.SetZ((clusterI->Z() - spdLength/2)*kconv); + rnew.SetQ(1.); + rnew.SetdEdX(0.); + rnew.SetSigmaX2(kRMSx*kRMSx); + rnew.SetSigmaZ2(kRMSz*kRMSz); + rnew.fTracks[0]=track0; + rnew.fTracks[1]=track1; + rnew.fTracks[2]=track2; + iTS->AddRecPoint(rnew); + } // I clusters + + fMap->ClearMap(); + +} +//_____________________________________________________________________________ + +void AliITSClusterFinderSPDdubna::FindRawClusters(Int_t mod) +{ + // find raw clusters + Find1DClusters(); + GroupClusters(); + TracksInCluster(); + GetRecPoints(); + +} + + + diff --git a/ITS/AliITSClusterFinderSPDdubna.h b/ITS/AliITSClusterFinderSPDdubna.h new file mode 100644 index 00000000000..906c84adb26 --- /dev/null +++ b/ITS/AliITSClusterFinderSPDdubna.h @@ -0,0 +1,69 @@ +#ifndef ALIITSCLUSTERFINDERSPDDUBNA_H +#define ALIITSCLUSTERFINDERSPDDUBNA_H + +//////////////////////////////////////////////// +// ITS Cluster Finder Class // +//////////////////////////////////////////////// + +#include "AliITSClusterFinder.h" + +class AliITSMapA1; + +class AliITSClusterFinderSPDdubna : + public AliITSClusterFinder + +{ +public: + AliITSClusterFinderSPDdubna + (AliITSsegmentation *seg, TClonesArray *dig, TClonesArray *recp); + AliITSClusterFinderSPDdubna(); + virtual ~AliITSClusterFinderSPDdubna(); + + void SetDx(Float_t dx=1.) { + // set dx + fDx=dx; + } + void SetDz(Float_t dz=0.) { + // set dz + fDz=dz; + } + void SetNCells(Int_t minc=0) { + // set ncells + fMinNCells=minc; + } + + // Search for clusters + void FindRawClusters(Int_t mod=0); + void Find1DClusters(); + void GroupClusters(); + void TracksInCluster(); + void SelectClusters() { + // selects clusters + } + void GetRecPoints(); + +private: + + AliITSClusterFinderSPDdubna(const AliITSClusterFinderSPDdubna &source); // copy ctor + AliITSClusterFinderSPDdubna& operator=(const AliITSClusterFinderSPDdubna &source); + +private: + + TClonesArray *fClusters; // clusters + Int_t fNclusters; // num of clusters + Float_t fDz; // dz + Float_t fDx; // dx + + Int_t fMinNCells; // min num of cells in the cluster + + ClassDef(AliITSClusterFinderSPDdubna,1) // SPD clustering - Boris B. algo based + // on Piergiorgio's algo + }; +#endif + + + + + + + diff --git a/ITS/AliITSFindClustersDubna.C b/ITS/AliITSFindClustersDubna.C new file mode 100644 index 00000000000..4380cb99679 --- /dev/null +++ b/ITS/AliITSFindClustersDubna.C @@ -0,0 +1,101 @@ +Int_t AliITSFindClustersDubna() { + + printf("FindClusters\n"); + + TFile *in=TFile::Open("galice.root","UPDATE"); + if (!in->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 2;} + + in->ls(); + + if (!(gAlice=(AliRun*)in->Get("gAlice"))) { + cerr<<"gAlice have not been found on galice.root !\n"; + return 2; + } + + + gAlice->GetEvent(0); + + AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS"); + if (!ITS) { + cerr<<"ITSFindClusters.C : AliITS object not found on file\n"; + return 3; + } + Int_t ver = ITS->IsVersion(); + cerr<<"ITS version "<MakeTreeC(); +// Set the models for cluster finding + AliITSgeom *geom = ITS->GetITSgeom(); + + // SPD + AliITSDetType *iDetType=ITS->DetType(0); + AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel(); + TClonesArray *dig0 = ITS->DigitsAddress(0); + TClonesArray *recp0 = ITS->ClustersAddress(0); + AliITSClusterFinderSPDdubna *rec0=new AliITSClusterFinderSPDdubna(seg0,dig0,recp0); + ITS->SetReconstructionModel(0,rec0); + + // SDD + AliITSDetType *iDetType=ITS->DetType(1); + AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel(); + if (!seg1) seg1 = new AliITSsegmentationSDD(geom); + AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel(); + if (!res1) res1=new AliITSresponseSDD(); + Float_t baseline,noise; + res1->GetNoiseParam(noise,baseline); + Float_t noise_after_el = res1->GetNoiseAfterElectronics(); + Float_t thres = baseline; + thres += (4.*noise_after_el); // TB // (4.*noise_after_el); + printf("thres %f\n",thres); + res1->Print(); + TClonesArray *dig1 = ITS->DigitsAddress(1); + TClonesArray *recp1 = ITS->ClustersAddress(1); + AliITSClusterFinderSDD *rec1=new AliITSClusterFinderSDD(seg1,res1,dig1,recp1); + rec1->SetCutAmplitude((int)thres); + ITS->SetReconstructionModel(1,rec1); + + + // SSD + AliITSDetType *iDetType=ITS->DetType(2); + AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel(); + TClonesArray *dig2 = ITS->DigitsAddress(2); + AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2); + ITS->SetReconstructionModel(2,rec2); + // test + printf("SSD dimensions %f %f \n",seg2->Dx(),seg2->Dz()); + printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx()); + + + + if(!gAlice->TreeR()) gAlice->MakeTree("R"); + //make branch + ITS->MakeBranch("R"); + + TStopwatch timer; + + switch (ver) { + case 5: + cerr<<"Looking for clusters...\n"; + { + timer.Start(); + ITS->DigitsToRecPoints(0,0,"All"); + } + break; + default: + cerr<<"Invalid ITS version !\n"; + return 5; + } + + timer.Stop(); timer.Print(); + + delete rec0; + delete rec1; + delete rec2; + + + delete gAlice; gAlice=0; + + in->Close(); + + return 0; +} diff --git a/ITS/AliITSHits2DigitsDubna.C b/ITS/AliITSHits2DigitsDubna.C new file mode 100644 index 00000000000..a571a24ad01 --- /dev/null +++ b/ITS/AliITSHits2DigitsDubna.C @@ -0,0 +1,109 @@ +Int_t AliITSHits2DigitsDubna() +{ + + // Connect the Root Galice file containing Geometry, Kine and Hits + + const char * inFile = "galice.root"; + TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(inFile); + if (file) {file->Close(); delete file;} + printf("Hits2Digits\n"); + file = new TFile(inFile,"UPDATE"); + if (!file->IsOpen()) { + cerr<<"Can't open "<ls(); + + // Get AliRun object from file or return if not on file + if (gAlice) delete gAlice; + gAlice = (AliRun*)file->Get("gAlice"); + if (!gAlice) { + cerr<<"ITSHits2Digits.C : AliRun object not found on file\n"; + return 2; + } + + gAlice->GetEvent(0); + AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS"); + if (!ITS) { + cerr<<"ITSHits2Digits.C : AliITS object not found on file\n"; + return 3; + } + +// Set the simulation models for the three detector types + AliITSgeom *geom = ITS->GetITSgeom(); + + // SPD + AliITSDetType *iDetType=ITS->DetType(0); + AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel(); + //AliITSresponseSPDdubna *res0 = (AliITSresponseSPDdubna*)iDetType->GetResponseModel(); + AliITSresponseSPDdubna *res0 = new AliITSresponseSPDdubna(); + ITS->SetResponseModel(0,res0); + AliITSsimulationSPDdubna *sim0=new AliITSsimulationSPDdubna(seg0,res0); + ITS->SetSimulationModel(0,sim0); + // test + printf("SPD dimensions %f %f \n",seg0->Dx(),seg0->Dz()); + printf("SPD npixels %d %d \n",seg0->Npz(),seg0->Npx()); + printf("SPD pitches %d %d \n",seg0->Dpz(0),seg0->Dpx(0)); + // end test + + // SDD + // Set response parameters + // SDD compression param: 2 fDecrease, 2fTmin, 2fTmax or disable, 2 fTolerance + AliITSDetType *iDetType=ITS->DetType(1); + AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel(); + if (!res1) { + res1=new AliITSresponseSDD(); + ITS->SetResponseModel(1,res1); + } + Float_t baseline; + Float_t noise; + res1->GetNoiseParam(noise,baseline); + Float_t noise_after_el = res1->GetNoiseAfterElectronics(); + cout << "noise_after_el: " << noise_after_el << endl; + Float_t fCutAmp; + fCutAmp = baseline; + fCutAmp += (2.*noise_after_el); // noise + cout << "Cut amplitude: " << fCutAmp << endl; + Int_t cp[8]={0,0,fCutAmp,fCutAmp,0,0,0,0}; + res1->SetCompressParam(cp); + + res1->Print(); + AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel(); + if (!seg1) { + seg1 = new AliITSsegmentationSDD(geom,res1); + ITS->SetSegmentationModel(1,seg1); + } + AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1); + ITS->SetSimulationModel(1,sim1); + + // SSD + AliITSDetType *iDetType=ITS->DetType(2); + AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel(); + AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel(); + res2->SetSigmaSpread(3.,2.); + AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2); + ITS->SetSimulationModel(2,sim2); + + cerr<<"Digitizing ITS...\n"; + + if(!gAlice->TreeD()) gAlice->MakeTree("D"); + printf("TreeD %p\n",gAlice->TreeD()); + //make branch + ITS->MakeBranch("D"); + + TStopwatch timer; + timer.Start(); + ITS->HitsToDigits(0,0,-1," ","All"," "); + timer.Stop(); timer.Print(); + + delete sim0; + delete sim1; + delete sim2; + + + delete gAlice; gAlice=0; + file->Close(); + delete file; + return 0; +}; + diff --git a/ITS/AliITSresponseSPDdubna.cxx b/ITS/AliITSresponseSPDdubna.cxx new file mode 100644 index 00000000000..bc6ee5f4089 --- /dev/null +++ b/ITS/AliITSresponseSPDdubna.cxx @@ -0,0 +1,32 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +#include + +#include "AliITSresponseSPDdubna.h" + +//___________________________________________ +ClassImp(AliITSresponseSPDdubna) + +AliITSresponseSPDdubna::AliITSresponseSPDdubna() +{ + // constructor + SetDiffCoeff(); + SetNoiseParam(); + SetDataType(); + SetMinVal(); + +} + diff --git a/ITS/AliITSresponseSPDdubna.h b/ITS/AliITSresponseSPDdubna.h new file mode 100644 index 00000000000..274d63b370e --- /dev/null +++ b/ITS/AliITSresponseSPDdubna.h @@ -0,0 +1,67 @@ +#ifndef ALIITSRESPONSESPDDUBNA_H +#define ALIITSRESPONSESPDDUBNA_H + +#include "AliITSresponse.h" +#include + +//---------------------------------------------- +// +// ITS response class for SPD +// +class AliITSresponseSPDdubna : + public AliITSresponse { +public: + + AliITSresponseSPDdubna(); + virtual ~AliITSresponseSPDdubna() { + // destructror + } + // + // Configuration methods + // + virtual void SetDiffCoeff(Float_t p1=0.00433,Float_t dummy=0.) { + // Diffusion coefficient + fDiffCoeff=p1; + } + virtual void DiffCoeff(Float_t &diffc,Float_t &dummy) { + // Get diffusion coefficient + diffc= fDiffCoeff; + } + virtual void SetNoiseParam(Float_t n=200., Float_t b=0.) { + // set noise and baseline + fNoise=n; fBaseline=b; + } + virtual void GetNoiseParam(Float_t &n, Float_t &b) { + // get noise and baseline + n=fNoise; b=fBaseline; + } + virtual void SetMinVal(Int_t p1=2000) { + // Zero-suppression option threshold + fThreshold=p1; + } + virtual Int_t MinVal() { + // Get zero-suppression threshold + return fThreshold; + } + virtual void SetDataType(const char *data="simulated") { + // Type of data - real or simulated + fDataType=data; + } + virtual const char *DataType() const { + // Get data typer + return fDataType.Data(); + } + + ClassDef(AliITSresponseSPDdubna,1) // SPD response + + protected: + + Float_t fDiffCoeff; // Diffusion Coefficient + Float_t fNoise; // Noise value + Float_t fBaseline; // Baseline value + Int_t fThreshold; // Zero-Suppression threshold + + TString fDataType; // Type of data - real or simulated +}; + +#endif diff --git a/ITS/AliITSsimulationSPDdubna.cxx b/ITS/AliITSsimulationSPDdubna.cxx new file mode 100644 index 00000000000..dd25d84a6fa --- /dev/null +++ b/ITS/AliITSsimulationSPDdubna.cxx @@ -0,0 +1,625 @@ +#include +#include +#include +#include +#include +#include + + +#include "AliRun.h" +#include "AliITS.h" +#include "AliITShit.h" +#include "AliITSdigit.h" +#include "AliITSmodule.h" +#include "AliITSMapA2.h" +#include "AliITSsimulationSPDdubna.h" +#include "AliITSsegmentation.h" +#include "AliITSresponse.h" + + + + +ClassImp(AliITSsimulationSPDdubna) +//////////////////////////////////////////////////////////////////////// +// Version: 0 +// Written by Boris Batyunya +// December 20 1999 +// +// AliITSsimulationSPDdubna is the simulation of SPDs +//________________________________________________________________________ + + +AliITSsimulationSPDdubna::AliITSsimulationSPDdubna() +{ + // constructor + fResponse = 0; + fSegmentation = 0; + fMapA2=0; + fHis = 0; + fNoise=0.; + fBaseline=0.; + fNPixelsZ=0; + fNPixelsX=0; +} + + +//_____________________________________________________________________________ + +AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(AliITSsegmentation *seg, AliITSresponse *resp) { + // standard constructor + + fHis = 0; + fResponse = resp; + fSegmentation = seg; + + fResponse->GetNoiseParam(fNoise,fBaseline); + + fMapA2 = new AliITSMapA2(fSegmentation); + + // + + fNPixelsZ=fSegmentation->Npz(); + fNPixelsX=fSegmentation->Npx(); + +} + +//_____________________________________________________________________________ + +AliITSsimulationSPDdubna::~AliITSsimulationSPDdubna() { + // destructor + + delete fMapA2; + + if (fHis) { + fHis->Delete(); + delete fHis; + } +} + + +//__________________________________________________________________________ +AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(const AliITSsimulationSPDdubna &source){ + // Copy Constructor + if(&source == this) return; + this->fMapA2 = source.fMapA2; + this->fNoise = source.fNoise; + this->fBaseline = source.fBaseline; + this->fNPixelsX = source.fNPixelsX; + this->fNPixelsZ = source.fNPixelsZ; + this->fHis = source.fHis; + return; +} + +//_________________________________________________________________________ +AliITSsimulationSPDdubna& + AliITSsimulationSPDdubna::operator=(const AliITSsimulationSPDdubna &source) { + // Assignment operator + if(&source == this) return *this; + this->fMapA2 = source.fMapA2; + this->fNoise = source.fNoise; + this->fBaseline = source.fBaseline; + this->fNPixelsX = source.fNPixelsX; + this->fNPixelsZ = source.fNPixelsZ; + this->fHis = source.fHis; + return *this; + } +//_____________________________________________________________________________ + +void AliITSsimulationSPDdubna::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy) +{ + // digitize module + + const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons + // for 3.6 eV/pair + const Float_t kconv = 10000.; // cm -> microns + + Float_t spdLength = fSegmentation->Dz(); + Float_t spdWidth = fSegmentation->Dx(); + + Float_t difCoef, dum; + fResponse->DiffCoeff(difCoef,dum); + + Float_t zPix0 = 1e+6; + Float_t xPix0 = 1e+6; + Float_t yPrev = 1e+6; + + Float_t zPitch = fSegmentation->Dpz(0); + Float_t xPitch = fSegmentation->Dpx(0); + + TObjArray *fHits = mod->GetHits(); + Int_t nhits = fHits->GetEntriesFast(); + if (!nhits) return; + + //cout<<"len,wid,dy,nx,nz,pitchx,pitchz ="<Dy()<<","<At(hit); + Int_t layer = iHit->GetLayer(); + Float_t yPix0 = -73; + if(layer == 1) yPix0 = -77; + + if(iHit->StatusEntering()) idhit=hit; + Int_t itrack = iHit->GetTrack(); + Int_t dray = 0; + + if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE; + + // Int_t parent = iHit->GetParticle()->GetFirstMother(); + Int_t partcode = iHit->GetParticle()->GetPdgCode(); + +// partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+ +// 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda + + /* + Float_t px = iHit->GetPXL(); // the momenta at the + Float_t py = iHit->GetPYL(); // each GEANT step + Float_t pz = iHit->GetPZL(); + Float_t ptot = 1000*sqrt(px*px+py*py+pz*pz); + */ + + Float_t pmod = iHit->GetParticle()->P(); // total momentum at the + // vertex + pmod *= 1000; + + + if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e- + // at p < 6 MeV/c + + + // Get hit z and x(r*phi) cordinates for each module (detector) + // in local system. + + Float_t zPix = kconv*iHit->GetZL(); + Float_t xPix = kconv*iHit->GetXL(); + Float_t yPix = kconv*iHit->GetYL(); + + // Get track status + Int_t status = iHit->GetTrackStatus(); + //cout<<"hit,status,y ="< spdLength/2) { + //cout<<"!!!1 z outside ="<MinVal(); + + Int_t digits[3], tracks[3], hits[3],gi,j1; + Float_t charges[3]; + Float_t electronics; + Float_t signal,phys; + for(Int_t iz=0;izGaus(); + signal = (float)fMapA2->GetSignal(iz,ix); + signal += electronics; + gi =iz*fNPixelsX+ix; // global index + if (signal > threshold) { + digits[0]=iz; + digits[1]=ix; + digits[2]=1; + for(j1=0;j1<3;j1++){ + if (pList[gi]) { + //b.b. tracks[j1]=-3; + tracks[j1] = (Int_t)(*(pList[gi]+j1)); + hits[j1] = (Int_t)(*(pList[gi]+j1+6)); + }else { + tracks[j1]=-2; //noise + hits[j1] = -1; + } + charges[j1] = 0; + } + + if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) { + tracks[1] = -3; + hits[1] = -1; + tracks[2] = -3; + hits[2] = -1; + } + if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) { + tracks[1] = -3; + hits[1] = -1; + } + if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) { + tracks[2] = -3; + hits[2] = -1; + } + if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) { + tracks[2] = -3; + hits[2] = -1; + } + + phys=0; + aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges); + } + if(pList[gi]) delete [] pList[gi]; + } + } + delete [] pList; + +} + + +//____________________________________________ + +void AliITSsimulationSPDdubna::CreateHistograms() +{ + // create 1D histograms for tests + + printf("SPD - create histograms\n"); + + fHis=new TObjArray(fNPixelsZ); + for (Int_t i=0;iReset(); + } + +} + + + + + + + + + diff --git a/ITS/AliITSsimulationSPDdubna.h b/ITS/AliITSsimulationSPDdubna.h new file mode 100644 index 00000000000..276138756da --- /dev/null +++ b/ITS/AliITSsimulationSPDdubna.h @@ -0,0 +1,51 @@ +#ifndef ALIITSSIMULATIONSPDDUBNA_H +#define ALIITSSIMULATIONSPDDUBNA_H + +#include "AliITSsimulation.h" + +class AliITSMapA2; +class AliITSsegmentation; +class AliITSresponse; +class AliITSmodule; + +//------------------------------------------------------------------- + +class AliITSsimulationSPDdubna : public AliITSsimulation { + +public: + + AliITSsimulationSPDdubna(); + AliITSsimulationSPDdubna(AliITSsegmentation *seg, AliITSresponse *res); + ~AliITSsimulationSPDdubna(); + AliITSsimulationSPDdubna(const AliITSsimulationSPDdubna &source); // copy constructor + AliITSsimulationSPDdubna& operator=(const AliITSsimulationSPDdubna &source); // ass. operator + + void DigitiseModule(AliITSmodule *mod,Int_t module,Int_t dummy); + void ChargeToSignal(Float_t **pList); + void GetList(Int_t track, Int_t hit, Float_t **pList, Int_t *IndexRange); + + void CreateHistograms(); + void ResetHistograms(); + TObjArray* GetHistArray() { + // get hist array + return fHis; + } + +private: + + AliITSMapA2 *fMapA2; // MapA2 + Float_t fNoise; // Noise + Float_t fBaseline; // Baseline + Int_t fNPixelsX; // NPixelsX + Int_t fNPixelsZ; // NPixelsZ + + TObjArray *fHis; // just in case for histogramming + + ClassDef(AliITSsimulationSPDdubna,1) // Simulation of SPD clusters + +}; + +#endif + + + diff --git a/ITS/AliITStestDubna.C b/ITS/AliITStestDubna.C new file mode 100644 index 00000000000..f959733bec1 --- /dev/null +++ b/ITS/AliITStestDubna.C @@ -0,0 +1,32 @@ +Int_t AliITStestDubna() { + Int_t rc=0; + +//Test ITS simulation + gROOT->LoadMacro("$(ALICE_ROOT)/macros/grun.C"); + grun(); + + Int_t ver=gAlice->GetDetector("ITS")->IsVersion(); + delete gAlice; gAlice=0; + + if (ver!=5) { + cerr<<"Invalid ITS version: "<LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2DigitsDubna.C"); + if (rc=AliITSHits2DigitsDubna()) return rc; + + } + + printf("start reconstruction\n"); + +//Test ITS reconstruction + gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersDubna.C"); + if (rc=AliITSFindClustersDubna()) return rc; + + //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSgraphycs.C"); + //if (rc=AliITSgraphycs()) return rc; + + return rc; +} diff --git a/ITS/ITSDigitsToClustersDubna.C b/ITS/ITSDigitsToClustersDubna.C new file mode 100644 index 00000000000..a71aa8dcf09 --- /dev/null +++ b/ITS/ITSDigitsToClustersDubna.C @@ -0,0 +1,161 @@ +#include "iostream.h" + +void ITSDigitsToClustersDubna (Int_t evNumber1=0,Int_t evNumber2=0) +{ +///////////////////////////////////////////////////////////////////////// +// This macro is a small example of a ROOT macro +// illustrating how to read the output of GALICE +// and do some analysis. +// +///////////////////////////////////////////////////////////////////////// + +// Dynamically link some shared libs + + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } else { + delete gAlice; + gAlice=0; + } + + +// Connect the Root Galice file containing Geometry, Kine and Hits + + TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); + printf("file %p\n",file); + if (file) file->Close(); + file = new TFile("galice.root","UPDATE"); + file->ls(); + + printf ("I'm after Map \n"); + +// Get AliRun object from file or create it if not on file + + if (!gAlice) { + gAlice = (AliRun*)file->Get("gAlice"); + if (gAlice) printf("AliRun object found on file\n"); + if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); + } + printf ("I'm after gAlice \n"); + + AliITS *ITS = (AliITS*) gAlice->GetModule("ITS"); + if (!ITS) return; + + AliITSgeom *geom = ITS->GetITSgeom(); + + + // NOTE: if you foresee to have (in segmentation or response) parameter + // values other than the default ones, and these values are used not only in + // simulation but in cluster finder as well, please set them via your + // local Config.C - the streamer will take care of writing the correct + // info and you'll no longer be obliged to set them again for your cluster + // finder as it's done in this macro + + + + // Set the models for cluster finding + + // SPD + + + + ITS->MakeTreeC(); + Int_t nparticles=gAlice->GetEvent(0); + + AliITSDetType *iDetType=ITS->DetType(0); + AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel(); + TClonesArray *dig0 = ITS->DigitsAddress(0); + TClonesArray *recp0 = ITS->ClustersAddress(0); + AliITSClusterFinderSPD *rec0=new AliITSClusterFinderSPD(seg0,dig0,recp0); + ITS->SetReconstructionModel(0,rec0); + + // SDD + + AliITSDetType *iDetType=ITS->DetType(1); + AliITSgeom *geom = ITS->GetITSgeom(); + + AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel(); + if (!seg1) seg1 = new AliITSsegmentationSDD(geom); + AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel(); + if (!res1) res1=new AliITSresponseSDD(); + + Float_t baseline,noise; + res1->GetNoiseParam(noise,baseline); + Float_t noise_after_el = res1->GetNoiseAfterElectronics(); + Float_t thres = baseline; + thres += (4.*noise_after_el); // TB // (4.*noise_after_el); + printf("thres %f\n",thres); + res1->Print(); + + TClonesArray *dig1 = ITS->DigitsAddress(1); + TClonesArray *recp1 = ITS->ClustersAddress(1); + AliITSClusterFinderSDD *rec1=new AliITSClusterFinderSDD(seg1,res1,dig1,recp1); + rec1->SetCutAmplitude((int)thres); + ITS->SetReconstructionModel(1,rec1); + rec1->Print(); + + // SSD + + AliITSDetType *iDetType=ITS->DetType(2); + AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel(); + seg2->SetDetSize(72960.,40000.,303.); + TClonesArray *dig2 = ITS->DigitsAddress(2); + AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2); + ITS->SetReconstructionModel(2,rec2); + // test + //printf("SSD dimensions %f %f \n",seg2->Dx(),seg2->Dz()); + //printf("SSD nstrips %d %d \n",seg2->Npz(),seg2->Npx()); + + + +// +// Event Loop +// + + cout << "Looking for clusters...\n"; + + TStopwatch timer; + + if(!gAlice->TreeR()) gAlice->MakeTree("R"); + //make branch + ITS->MakeBranch("R"); + + for (int nev=evNumber1; nev<= evNumber2; nev++) { + if(nev>0) { + nparticles = gAlice->GetEvent(nev); + gAlice->SetEvent(nev); + if(!gAlice->TreeR()) gAlice-> MakeTree("R"); + ITS->MakeBranch("R"); + } + cout << "nev " <DigitsToRecPoints(nev,last_entry,"All"); + timer.Stop(); timer.Print(); + } // event loop + + delete rec0; + delete rec1; + delete rec2; + + file->Close(); +} + + + + + + + + + + + + + + diff --git a/ITS/ITSHitsToDigitsDubna.C b/ITS/ITSHitsToDigitsDubna.C new file mode 100644 index 00000000000..1146de9d10e --- /dev/null +++ b/ITS/ITSHitsToDigitsDubna.C @@ -0,0 +1,177 @@ +#include "iostream.h" + +void ITSHitsToDigitsDubna (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nsignal =25, Int_t size=-1) +{ +///////////////////////////////////////////////////////////////////////// +// This macro is a small example of a ROOT macro +// illustrating how to read the output of GALICE +// and do some analysis. +// +///////////////////////////////////////////////////////////////////////// + +// Dynamically link some shared libs + + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } else { + delete gAlice; + gAlice=0; + } + + +// Connect the Root Galice file containing Geometry, Kine and Hits + + TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); + printf("file %p\n",file); + if (file) file->Close(); + if (!file) file = new TFile("galice.root","UPDATE"); + file->ls(); + + printf ("I'm after Map \n"); + +// Get AliRun object from file or create it if not on file + + if (!gAlice) { + gAlice = (AliRun*)file->Get("gAlice"); + if (gAlice) printf("AliRun object found on file\n"); + if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); + } + printf ("I'm after gAlice \n"); + + AliITS *ITS = (AliITS*) gAlice->GetModule("ITS"); + if (!ITS) return; + + + // Set the simulation models + + AliITSgeom *geom = ITS->GetITSgeom(); + + // SDD + // SDD compression param: 2 fDecrease, 2fTmin, 2fTmax or disable, 2 fTolerance + + AliITSDetType *iDetType=ITS->DetType(1); + AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel(); + if (!res1) { + res1=new AliITSresponseSDD(); + ITS->SetResponseModel(1,res1); + } + + //res1->SetChargeLoss(0.); + Float_t baseline; + Float_t noise; + res1->GetNoiseParam(noise,baseline); + Float_t noise_after_el = res1->GetNoiseAfterElectronics(); + cout << "noise_after_el: " << noise_after_el << endl; + Float_t fCutAmp; + fCutAmp = baseline; + fCutAmp += (2.*noise_after_el); // noise + cout << "Cut amplitude: " << fCutAmp << endl; + Int_t cp[8]={0,0,fCutAmp,fCutAmp,0,0,0,0}; + res1->SetCompressParam(cp); + // res1->SetElectronics(2); // 1 = Pascal, 2 = OLA + + res1->Print(); + + //cout << "SDD segmentation" << endl; + + AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel(); + if (!seg1) { + seg1 = new AliITSsegmentationSDD(geom,res1); + ITS->SetSegmentationModel(1,seg1); + } + seg1->Print(); + + //cout << "SDD segmentation" << endl; + AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1); + ITS->SetSimulationModel(1,sim1); + sim1->Print(); + + + + // SPD + + AliITSDetType *iDetType=ITS->DetType(0); + AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel(); + //AliITSresponseSPDdubna *res0 = (AliITSresponseSPDdubna*)iDetType->GetResponseModel(); + AliITSresponseSPDdubna *res0 = new AliITSresponseSPDdubna(); + ITS->SetResponseModel(0,res0); + AliITSsimulationSPDdubna *sim0=new AliITSsimulationSPDdubna(seg0,res0); + ITS->SetSimulationModel(0,sim0); + // test + //printf("SPD dimensions %f %f \n",seg0->Dx(),seg0->Dz()); + //printf("SPD npixels %d %d \n",seg0->Npz(),seg0->Npx()); + //printf("SPD pitches %d %d \n",seg0->Dpz(0),seg0->Dpx(0)); + // end test + + + // SSD + + AliITSDetType *iDetType=ITS->DetType(2); + AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel(); + seg2->SetDetSize(72960.,40000.,303.); + AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel(); + res2->SetSigmaSpread(3.,2.); + AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2); + ITS->SetSimulationModel(2,sim2); + + +// +// Event Loop +// + + + // create the TreeD + + Int_t nparticles=gAlice->GetEvent(0); + printf("Create TreeD \n"); + if(!gAlice->TreeD()) gAlice->MakeTree("D"); + //make branch + ITS->MakeBranch("D"); + + Int_t nbgr_ev=0; + + + cout<<"Digitizing ITS...\n"; + TStopwatch timer; + + for (Int_t nev=evNumber1; nev<= evNumber2; nev++) { + cout << "nev " <0) { + nparticles = gAlice->GetEvent(nev); + gAlice->SetEvent(nev); + if(!gAlice->TreeD()) gAlice-> MakeTree("D"); + ITS->MakeBranch("D"); + } + cout << "nparticles " <HitsToDigits(nev,nbgr_ev,size," ","All"," "); + timer.Stop(); timer.Print(); + } // event loop + + delete sim0; + delete sim1; + delete sim2; + + + file->Close(); +} + + + + + + + + + + + + + + diff --git a/ITS/ITSreadClustTestSPDDubna.C b/ITS/ITSreadClustTestSPDDubna.C new file mode 100644 index 00000000000..e79ded66424 --- /dev/null +++ b/ITS/ITSreadClustTestSPDDubna.C @@ -0,0 +1,111 @@ +#include "iostream.h" + +void ITSreadClustTestSPDDubna (Int_t evNumber1=0,Int_t evNumber2=0) +{ +///////////////////////////////////////////////////////////////////////// +// This macro is a small example of a ROOT macro +// illustrating how to read the output of GALICE +// and do some analysis. +// +///////////////////////////////////////////////////////////////////////// + +// Dynamically link some shared libs + + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } + // Anyway, this macro needs to read a gAlice file, so it + // clears the gAlice object if there is already one in memory... + else { + if(gAlice){ + delete gAlice; + gAlice = 0; + } + } + +// Connect the Root Galice file containing Geometry, Kine and Hits + + TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); + if (!file) file = new TFile("galice.root"); + file->ls(); + +// Get AliRun object from file or create it if not on file + + if (!gAlice) { + gAlice = (AliRun*)file->Get("gAlice"); + if (gAlice) printf("AliRun object found on file\n"); + if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); + } + +// +// Loop over events +// + Int_t Nh=0; + Int_t Nh1=0; + for (int nev=0; nev<= evNumber2; nev++) { + Int_t nparticles = gAlice->GetEvent(nev); + cout << "nev " << nev <TreeH(); + Int_t ntracks = TH->GetEntries(); + cout<<"ntracks "<GetModule("ITS"); + if(!ITS) return; + TClonesArray *Particles = gAlice->Particles(); + + ITS->GetTreeC(nev); + TTree *TC=ITS->TreeC(); + Int_t nent=TC->GetEntries(); + printf("Found %d entries in the tree (must be one per module per event!)\n",nent); + + for (Int_t idet=0;idet<3;idet++) { + + TClonesArray *ITSclu = ITS->ClustersAddress(idet); + + if (idet != 0) continue; + for (Int_t mod=0; modResetClusters(); + TC->GetEvent(mod); + Int_t ncl = ITSclu->GetEntries(); + if (ncl) printf("Found %d clusters for module %d in det type %d \n",ncl,mod,idet); + + if (!ncl) continue; + + for (Int_t icl=0;iclUncheckedAt(icl); + printf("%d %d %f %f %f\n",ITSclust->NclZ(),ITSclust->NclX(),ITSclust->Q(),ITSclust->X(),ITSclust->Z()); + } + //for (Int_t icl=0;iclUncheckedAt(icl); + //printf("%d %d %f %f %f\n",ITSclust->Anodes(),ITSclust->Samples(),ITSclust->Q(),ITSclust->X(),ITSclust->Z()); + //} + } + } + + } // event loop + + cout<<"END test for clusters "<Close(); +} + + + + + + + + + + diff --git a/ITS/SPDclusterTestDubna.C b/ITS/SPDclusterTestDubna.C new file mode 100644 index 00000000000..fefcd0254db --- /dev/null +++ b/ITS/SPDclusterTestDubna.C @@ -0,0 +1,565 @@ +#include "iostream.h" + +void SPDclusterTestDubna (Int_t evNumber1=0,Int_t evNumber2=0) +{ +///////////////////////////////////////////////////////////////////////// +// This macro is a small example of a ROOT macro +// illustrating how to read the output of GALICE +// and do some analysis. +// +///////////////////////////////////////////////////////////////////////// + +// Dynamically link some shared libs + + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } else { + delete gAlice; + gAlice=0; + } + +// Connect the Root Galice file containing Geometry, Kine and Hits + + TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); + if (!file) file = new TFile("galice.root"); + file->ls(); + +// Get AliRun object from file or create it if not on file + + if (!gAlice) { + gAlice = (AliRun*)file->Get("gAlice"); + if (gAlice) printf("AliRun object found on file\n"); + if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); + } + + // ------------ Cluster and point analysis histogramms ------------ + + TH1F *Nxpix1 = new TH1F("Nxpix1","Cluster size in x(r*phi) direction for layer 1",20,0.,20.); + TH1F *Nxpix2 = new TH1F("Nxpix2","Cluster size in x(r*phi) direction for layer 2",20,0.,20.); + TH1F *Nzpix1 = new TH1F("Nzpix1","Cluster size in z direction for layer 1",15,0.,15.); + TH1F *Nzpix2 = new TH1F("Nzpix2","Cluster size in z direction for layer 2",15,0.,15.); + TH1F *Xpix1 = new TH1F("Xpix1","Local x coordinate (mm) for layer 1",20,-2.,18.); + TH1F *Xpix2 = new TH1F("Xpix2","Local x coordinate (mm) for layer 2",20,-2.,18.); + TH1F *Zpix1 = new TH1F("Zpix1","Local z coordinate (mm) for layer 1",90,-2.,88.); + TH1F *Zpix2 = new TH1F("Zpix2","Lolac z coordinate (mm) for layer 2",90,-2.,88.); + + TH1F *Xres1 = new TH1F("Xres1","Xrec and Xgen difference (micr) for layers 1",100,-200.,200.); + TH1F *Xres2 = new TH1F("Xres2","Xrec and Xgen difference (micr) for layers 2",100,-200.,200.); + TH1F *Zres1 = new TH1F("Zres1","Zrec and Zgen difference (micr) for layers 1",100,-800.,800.); + TH1F *Zres2 = new TH1F("Zres2","Zrec and Zgen difference (micr) for layers 2",100,-800.,800.); + + TH1F *Ptot1 = new TH1F("Ptot1","Total momentum (GeV/C) for layers 1",100,0.,5.); + TH1F *Pz1 = new TH1F("Pz1","Pz (GeV/C) for layers 1",100,-5.,5.); + TH1F *Theta1 = new TH1F("Theta1","Theta angle (rad) for layers 1",100,0.,4.); + TH1F *Y1 = new TH1F("Y1","Rapidity for layers 1",100,-4.,4.); + TH1F *Eta1 = new TH1F("Eta1","PseudoRapidity for layers 1",100,-4.,4.); + TH1F *Y1Den = new TH1F("Y1Den","Rapidity for layers 1",100,-0.5,0.5); + TH1F *Eta1Den = new TH1F("Eta1Den","PseudoRapidity for layers 1",100,-0.5,0.5); + TH1F *Y1DenA = new TH1F("Y1DenA","Rapidity for layers 1",100,-0.5,0.5); + TH1F *Eta1DenA = new TH1F("Eta1DenA","PseudoRapidity for layers 1",100,-0.5,0.5); + TH1F *Phi1 = new TH1F("Phi1","Phi angle (rad) for layers 1",100,0.,7.); + TH1F *Ptot2 = new TH1F("Ptot2","Total momentum (GeV/C) for layers 2",100,0.,5.); + TH1F *Pz2 = new TH1F("Pz2","Pz (GeV/C) for layers 2",100,-5.,5.); + TH1F *Theta2 = new TH1F("Theta2","Theta angle (rad) for layers 2",100,0.,4.); + TH1F *Y2 = new TH1F("Y2","Rapidity for layers 2",100,-4.,4.); + TH1F *Eta2 = new TH1F("Eta2","PseudoRapidity for layers 2",100,-4.,4.); + TH1F *Y2Den = new TH1F("Y2Den","Rapidity for layers 2",100,-0.5,0.5); + TH1F *Eta2Den = new TH1F("Eta2Den","PseudoRapidity for layers 2",100,-0.5,0.5); + TH1F *Y2DenA = new TH1F("Y2DenA","Rapidity for layers 2",100,-0.5,0.5); + TH1F *Eta2DenA = new TH1F("Eta2DenA","PseudoRapidity for layers 2",100,-0.5,0.5); + TH1F *Phi2 = new TH1F("Phi2","Phi angle (rad) for layers 2",100,0.,7.); + + // -------------- Create ntuples -------------------- + // ntuple structures: + + struct { + Int_t lay; + Int_t nx; + Int_t nz; + Int_t hitprim; + Int_t partcode; + Int_t ntrover; + Float_t dx; + Float_t dz; + Float_t pmod; + } ntuple_st; + + struct { + Int_t lay; + Int_t lad; + Int_t det; + Int_t nx; + Int_t nz; + Int_t ntrover; + Int_t noverlaps; + Int_t noverprim; + Float_t qcl; + Float_t dx; + Float_t dz; + Float_t x; + Float_t z; + } ntuple1_st; + + struct { + // Int_t lay; + Int_t lay; + Int_t nx; + Int_t nz; + Float_t x; + Float_t z; + Float_t qcl; + } ntuple2_st; + + ntuple = new TTree("ntuple","Demo ntuple"); + ntuple->Branch("lay",&ntuple_st.lay,"lay/I"); + ntuple->Branch("nx",&ntuple_st.nx,"nx/I"); + ntuple->Branch("nz",&ntuple_st.nz,"nz/I"); + ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I"); + ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I"); + ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I"); + ntuple->Branch("dx",&ntuple_st.dx,"dx/F"); + ntuple->Branch("dz",&ntuple_st.dz,"dz/F"); + ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F"); + + ntuple1 = new TTree("ntuple1","Demo ntuple1"); + ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I"); + ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I"); + ntuple1->Branch("det",&ntuple1_st.det,"det/I"); + ntuple1->Branch("nx",&ntuple1_st.nx,"nx/I"); + ntuple1->Branch("nz",&ntuple1_st.nz,"nz/I"); + ntuple1->Branch("qcl",&ntuple1_st.qcl,"qcl/F"); + ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I"); + ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I"); + ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I"); + ntuple1->Branch("x",&ntuple1_st.x,"x/F"); + ntuple1->Branch("z",&ntuple1_st.z,"z/F"); + ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F"); + ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F"); + + + ntuple2 = new TTree("ntuple2","Demo ntuple2"); + // ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I"); + ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I"); + ntuple2->Branch("x",&ntuple2_st.x,"x/F"); + ntuple2->Branch("z",&ntuple2_st.z,"z/F"); + ntuple2->Branch("nx",&ntuple2_st.nx,"nx/I"); + ntuple2->Branch("nz",&ntuple2_st.nz,"nz/I"); + ntuple2->Branch("qcl",&ntuple2_st.qcl,"qcl/F"); + +// ------------------------------------------------------------------------ +// +// Loop over events +// + for (int nev=0; nev<= evNumber2; nev++) { + Int_t nparticles = gAlice->GetEvent(nev); + cout << "nev " << nev <TreeH(); + Int_t ntracks = TH->GetEntries(); + cout<<"ntracks "<GetModule("ITS"); + TClonesArray *Particles = gAlice->Particles(); + + if (ITS) { + // fill modules with sorted by module hits + Int_t nmodules; + ITS->InitModules(-1,nmodules); + // ITS->FillModules(nev,-1,evNumber2,nmodules," "," "); + ITS->FillModules(nev,evNumber2,nmodules," "," "); + //get pointer to modules array + TObjArray *ITSmodules = ITS->GetModules(); + AliITShit *itsHit; + + // get the Tree for clusters + ITS->GetTreeC(nev); + TTree *TC=ITS->TreeC(); + Int_t nent=TC->GetEntries(); + printf("Found %d entries in the tree (must be one per module per event!)\n",nent); + Int_t lay, lad, det; + AliITSgeom *geom = ITS->GetITSgeom(); + + for (Int_t idettype=0;idettype<3;idettype++) { + + TClonesArray *ITSclusters = ITS->ClustersAddress(idettype); + //printf ("ITSclusters %p \n",ITSclusters); + + if (idettype != 0) continue; + + Float_t occup1 = 0; + Float_t occup2 = 0; + + // Module loop + for (Int_t mod=0; modAt(mod); + geom->GetModuleId(mod,lay,lad,det); + + Int_t nhits = itsModule->GetNhits(); + //if(nhits) printf("module nhits %d %d\n",mod,nhits); + if(!nhits) continue; + + ITS->ResetClusters(); + TC->GetEvent(mod); + Int_t nclust = ITSclusters->GetEntries(); + if (!nclust) continue; + + // cluster/hit loops + //cout<<"mod,lay,nclust,nhits ="<UncheckedAt(clu); + + Int_t noverlaps = 0; + Int_t noverprim = 0; + + Int_t clustersizex = itsclu->NclX(); + Int_t clustersizez = itsclu->NclZ(); + // Int_t xstart = itsclu->XStart(); + // Int_t xstop = itsclu->XStop(); + Int_t xstart = itsclu->XStartf(); + Int_t xstop = itsclu->XStopf(); + Float_t fxstart = xstart*50; + Float_t fxstop = (xstop+1)*50; + Float_t zstart = itsclu->ZStart(); + Float_t zstop = itsclu->ZStop(); + Int_t zend = itsclu->Zend(); + Int_t ntrover = itsclu->NTracks(); + Float_t clusterx = itsclu->X(); + Float_t clusterz = itsclu->Z(); + Float_t clusterQ = itsclu->Q(); + + if(lay == 1) occup1 += clusterQ; + if(lay == 2) occup2 += clusterQ; + + ntuple2_st.lay = lay; + ntuple2_st.x = clusterx/1000.; + ntuple2_st.z = clusterz/1000.; + ntuple2_st.nx = clustersizex; + ntuple2_st.nz = clustersizez; + ntuple2_st.qcl = clusterQ; + + ntuple2->Fill(); + + Int_t icl = 0; + Float_t dxprimlast = 10.e+6; + Float_t dzprimlast = 10.e+6; + + Float_t SPDlength = 83600; + Float_t SPDwidth = 12800; + Float_t xhit0 = 1e+5; + Float_t zhit0 = 1e+5; + + for (Int_t hit=0;hitGetHit(hit); + + Int_t hitlayer = itsHit->GetLayer(); + Int_t hitladder= itsHit->GetLadder(); + Int_t hitdet= itsHit->GetDetector(); + Float_t dEn = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV + Int_t track = itsHit->GetTrack(); + Int_t dray = 0; + Int_t hitstat = itsHit->GetTrackStatus(); + + Float_t zhit = 10000*itsHit->GetZL(); + Float_t xhit = 10000*itsHit->GetXL(); + + Float_t pxsimL = itsHit->GetPXL(); // the momenta at GEANT points + Float_t pysimL = itsHit->GetPYL(); + Float_t pzsimL = itsHit->GetPZL(); + Float_t psimL = TMath::Sqrt(pxsimL*pxsimL+pysimL*pysimL+pzsimL*pzsimL); + + // Check boundaries + if(zhit > SPDlength/2) { + //cout<<"!!! z outside ="<GetPXG(); // the momenta at this GEANT point + Float_t pysim = itsHit->GetPYG(); + Float_t pzsim = itsHit->GetPZG(); + Float_t psim = TMath::Sqrt(pxsim*pxsim+pysim*pysim+pzsim*pzsim); + + Int_t hitprim = 0; + + if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e- + // at p < 6 MeV/c + + if(dray == 0) noverlaps = noverlaps + 1; // overlapps for all hits but + // not for delta ray which + // also went out from the + // detector and returned + // again + + if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery + // particles + + if(hitprim > 0) noverprim = noverprim + 1; + + if(hitprim > 0) { + dxprimlast = xdif; + dzprimlast = zdif; + } + + // fill ntuple + + ntuple_st.lay = hitlayer; + ntuple_st.nx = clustersizex; + ntuple_st.nz = clustersizez; + ntuple_st.hitprim = hitprim; + ntuple_st.partcode = partcode; + ntuple_st.ntrover = ntrover; + ntuple_st.dx = xdif; + ntuple_st.dz = zdif; + ntuple_st.pmod = pmod; + + ntuple->Fill(); + + if(hitlayer == 1) { + Y1DenA->Fill(y); + Eta1DenA->Fill(eta); + } + if(hitlayer == 2) { + Y2DenA->Fill(y); + Eta2DenA->Fill(eta); + } + + + + if(hitprim > 0) { // for primary particles + + if(hitlayer == 1) { + Xres1->Fill(xdif); + Zres1->Fill(zdif); + Ptot1->Fill(pmod/1000.); + Pz1->Fill(pz); + Theta1->Fill(theta); + Y1->Fill(y); + Eta1->Fill(eta); + Y1Den->Fill(y); + Eta1Den->Fill(eta); + Phi1->Fill(phi); + } + if(hitlayer == 2) { + Xres2->Fill(xdif); + Zres2->Fill(zdif); + Ptot2->Fill(pmod/1000.); + Pz2->Fill(pz); + Theta2->Fill(theta); + Y2->Fill(y); + Eta2->Fill(eta); + Y2Den->Fill(y); + Eta2Den->Fill(eta); + Phi2->Fill(phi); + } + } // primery particles + + } // end of cluster region + } // end of hit loop + + if(icl == 1) { + + // fill ntuple1 + + // ntuple1->Fill(clusterlayer,clustersizex,clustersizez,noverlaps,\ +noverprim,dx,dz); + + if(noverlaps == 0) noverlaps = 1; // cluster contains one or more + // delta rays only + + ntuple1_st.lay = lay; + ntuple1_st.lad = lad; + ntuple1_st.det = det; + ntuple1_st.x = clusterx*1000.; + ntuple1_st.z = clusterz*1000.; + ntuple1_st.nx = clustersizex; + ntuple1_st.nz = clustersizez; + ntuple1_st.qcl = clusterQ; + ntuple1_st.ntrover = ntrover; + ntuple1_st.noverlaps = noverlaps; + ntuple1_st.noverprim = noverprim; + ntuple1_st.dx = dxprimlast; + ntuple1_st.dz = dzprimlast; + + ntuple1->Fill(); + + } // icl = 1 + } // cluster loop + } // module loop + + cout<<" Occupancy for layer-1 ="<Write(); + ntuple1->Write(); + ntuple2->Write(); + + Nxpix1->Write(); + Nzpix1->Write(); + Nxpix2->Write(); + Nzpix2->Write(); + + Xpix1->Write(); + Zpix1->Write(); + Xpix2->Write(); + Zpix2->Write(); + + Xres1->Write(); + Zres1->Write(); + Xres2->Write(); + Zres2->Write(); + + Ptot1->Write(); + Pz1->Write(); + Theta1->Write(); + Y1->Write(); + Eta1->Write(); + Y1Den->Write(); + Eta1Den->Write(); + Y1DenA->Write(); + Eta1DenA->Write(); + Phi1->Write(); + + Ptot2->Write(); + Pz2->Write(); + Theta2->Write(); + Y2->Write(); + Eta2->Write(); + Y2Den->Write(); + Eta2Den->Write(); + Y2DenA->Write(); + Eta2DenA->Write(); + Phi2->Write(); + + fhistos.Close(); + cout<<"!!! Histogramms and ntuples were written"<Divide(2,2); + c1->cd(1); + gPad->SetFillColor(33); + Xres1->SetFillColor(42); + Xres1->Draw(); + c1->cd(2); + gPad->SetFillColor(33); + Zres1->SetFillColor(46); + Zres1->Draw(); + c1->cd(3); + gPad->SetFillColor(33); + Xres2->SetFillColor(42); + Xres2->Draw(); + c1->cd(4); + gPad->SetFillColor(33); + Zres2->SetFillColor(46); + Zres2->Draw(); + + cout<<"END test for clusters and hits "<Close(); +} + + + -- 2.43.0