--- /dev/null
+/**************************************************************************
+ * 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 <iostream.h>
+#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<fNofPixels;k++) {
+
+ Int_t mmax = 10; // a size of the window for the cluster finding
+
+ for(it=0;it<fMaxNofSamples;it++) {
+
+ Int_t lclx = 0;
+ Int_t xstart = 0;
+ Int_t xstop = 0;
+ Int_t id = 0;
+ Int_t ilcl =0;
+
+ for(m=0;m<mmax;m++) { // find the cluster inside the window
+ id = it+m;
+ if(id >= 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; i<nofClusters; i++) label[i] = 0;
+ for(i=0; i<nofClusters; i++) {
+ if(label[i] != 0) continue;
+ for(j=i+1; j<nofClusters; j++) {
+ if(label[j] != 0) continue;
+ clusterI = (AliITSRawClusterSPD*) fClusters->At(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"<<endl;
+
+ /*
+ cout << "clusters " << i << "," << j << " before grouping" << endl;
+ clusterI->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<<endl;
+ */
+
+ delete [] label;
+
+ return;
+
+
+}
+//_____________________________________________________________________________
+
+void AliITSClusterFinderSPDdubna::TracksInCluster()
+{
+
+ // Find tracks creating one cluster
+
+ // get number of clusters for this module
+ Int_t nofClusters = fClusters->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; i<nofClusters; i++) {
+ ii = 0;
+ memset(cltracks,-1,sizeof(int)*trmax);
+ tr0=tr1=tr2=-1;
+
+ AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(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; iz<nclz; iz++) {
+ jz = zstart + iz;
+ for(ix=0; ix<nclx; ix++) {
+ jx = xstart + ix;
+ ind = fMap->GetHitIndex(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<trmax; is++) {
+ if(cltracks[is]<0) continue;
+ for(js=is+1; js<trmax; js++) {
+ if(cltracks[js]<0) continue;
+ if(cltracks[js]==cltracks[is]) cltracks[js]=-5;
+ }
+ }
+
+ ntr = 0;
+ for(ie=0; ie<trmax; ie++) {
+ if(cltracks[ie] >= 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; i<nofClusters; i++) {
+
+ AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(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();
+
+}
+
+
+
--- /dev/null
+#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
+
+
+
+
+
+
+
--- /dev/null
+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 "<<ver<<" has been found !\n";
+
+ ITS->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;
+}
--- /dev/null
+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 "<<inFile<<" !\n";
+ return 1;
+ }
+ file->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;
+};
+
--- /dev/null
+/**************************************************************************
+ * 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 <TMath.h>
+
+#include "AliITSresponseSPDdubna.h"
+
+//___________________________________________
+ClassImp(AliITSresponseSPDdubna)
+
+AliITSresponseSPDdubna::AliITSresponseSPDdubna()
+{
+ // constructor
+ SetDiffCoeff();
+ SetNoiseParam();
+ SetDataType();
+ SetMinVal();
+
+}
+
--- /dev/null
+#ifndef ALIITSRESPONSESPDDUBNA_H
+#define ALIITSRESPONSESPDDUBNA_H
+
+#include "AliITSresponse.h"
+#include <TString.h>
+
+//----------------------------------------------
+//
+// 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
--- /dev/null
+#include <iostream.h>
+#include <TRandom.h>
+#include <TH1.h>
+#include <TMath.h>
+#include <TString.h>
+#include <TParticle.h>
+
+
+#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 ="<<spdLength<<","<<spdWidth<<","<<fSegmentation->Dy()<<","<<fNPixelsX<<","<<fNPixelsZ<<","<<xPitch<<","<<zPitch<<endl;
+ // Array of pointers to the label-signal list
+
+ Int_t maxNDigits = fNPixelsX*fNPixelsZ + fNPixelsX ;;
+ Float_t **pList = new Float_t* [maxNDigits];
+ memset(pList,0,sizeof(Float_t*)*maxNDigits);
+ Int_t indexRange[4] = {0,0,0,0};
+
+ // Fill detector maps with GEANT hits
+ // loop over hits in the module
+ static Bool_t first;
+ Int_t lasttrack=-2;
+ Int_t hit, iZi, jz, jx;
+ //cout<<"SPD: module,nhits ="<<module<<","<<nhits<<endl;
+ Int_t idhit=-1;
+ for (hit=0;hit<nhits;hit++) {
+ AliITShit *iHit = (AliITShit*) fHits->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 ="<<hit<<","<<status<<","<<yPix<<endl;
+
+ // Check boundaries
+ if(zPix > spdLength/2) {
+ //cout<<"!!!1 z outside ="<<zPix<<endl;
+ zPix = spdLength/2 - 10;
+ //cout<<"!!!2 z outside ="<<zPix<<endl;
+ }
+ if(zPix < 0 && zPix < -spdLength/2) {
+ //cout<<"!!!1 z outside ="<<zPix<<endl;
+ zPix = -spdLength/2 + 10;
+ //cout<<"!!!2 z outside ="<<zPix<<endl;
+ }
+ if(xPix > spdWidth/2) {
+ //cout<<"!!!1 x outside ="<<xPix<<endl;
+ xPix = spdWidth/2 - 10;
+ //cout<<"!!!2 x outside ="<<xPix<<endl;
+ }
+ if(xPix < 0 && xPix < -spdWidth/2) {
+ //cout<<"!!!1 x outside ="<<xPix<<endl;
+ xPix = -spdWidth/2 + 10;
+ //cout<<"!!!2 x outside ="<<xPix<<endl;
+ }
+ Int_t trdown = 0;
+
+ // enter Si or after event in Si
+ if (status == 66 ) {
+ zPix0 = zPix;
+ xPix0 = xPix;
+ yPrev = yPix;
+ }
+
+ Float_t depEnergy = iHit->GetIonization();
+ // skip if the input point to Si
+
+ if(depEnergy <= 0.) continue;
+
+ // if track returns to the opposite direction:
+ if (yPix < yPrev) {
+ trdown = 1;
+ }
+
+
+ // take into account the holes diffusion inside the Silicon
+ // the straight line between the entrance and exit points in Si is
+ // divided into the several steps; the diffusion is considered
+ // for each end point of step and charge
+ // is distributed between the pixels through the diffusion.
+
+
+ // ---------- the diffusion in Z (beam) direction -------
+
+ Float_t charge = depEnergy*kEnToEl; // charge in e-
+ Float_t drPath = 0.;
+ Float_t tang = 0.;
+ Float_t sigmaDif = 0.;
+ Float_t zdif = zPix - zPix0;
+ Float_t xdif = xPix - xPix0;
+ Float_t ydif = TMath::Abs(yPix - yPrev);
+ Float_t ydif0 = TMath::Abs(yPrev - yPix0);
+
+ if(ydif < 1) continue; // ydif is not zero
+
+ Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
+
+ Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
+ Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1;
+
+ // number of the steps along the track:
+ Int_t nsteps = ndZ;
+ if(ndX > ndZ) nsteps = ndX;
+ if(nsteps < 6) nsteps = 6; // minimum number of the steps
+
+ if (projDif < 5 ) {
+ drPath = (yPix-yPix0)*1.e-4;
+ drPath = TMath::Abs(drPath); // drift path in cm
+ sigmaDif = difCoef*sqrt(drPath); // sigma diffusion in cm
+ sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
+ nsteps = 1;
+ }
+
+ if(projDif > 5) tang = ydif/projDif;
+ Float_t dCharge = charge/nsteps; // charge in e- for one step
+ Float_t dZ = zdif/nsteps;
+ Float_t dX = xdif/nsteps;
+
+ for (iZi = 1;iZi <= nsteps;iZi++) {
+ Float_t dZn = iZi*dZ;
+ Float_t dXn = iZi*dX;
+ Float_t zPixn = zPix0 + dZn;
+ Float_t xPixn = xPix0 + dXn;
+
+ if(projDif >= 5) {
+ Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
+ drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm
+ if(trdown == 0) {
+ drPath = TMath::Abs(drPath) + ydif0*1.e-4;
+ }
+ if(trdown == 1) {
+ drPath = ydif0*1.e-4 - TMath::Abs(drPath);
+ drPath = TMath::Abs(drPath);
+ }
+ sigmaDif = difCoef*sqrt(drPath);
+ sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
+ }
+
+ zPixn = (zPixn + spdLength/2.);
+ xPixn = (xPixn + spdWidth/2.);
+ Int_t nZpix, nXpix;
+ fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
+ zPitch = fSegmentation->Dpz(nZpix);
+ fSegmentation->GetPadTxz(xPixn,zPixn);
+ // set the window for the integration
+ Int_t jzmin = 1;
+ Int_t jzmax = 3;
+ if(nZpix == 1) jzmin =2;
+ if(nZpix == fNPixelsZ) jzmax = 2;
+
+ Int_t jxmin = 1;
+ Int_t jxmax = 3;
+ if(nXpix == 1) jxmin =2;
+ if(nXpix == fNPixelsX) jxmax = 2;
+
+ Float_t zpix = nZpix;
+ Float_t dZright = zPitch*(zpix - zPixn);
+ Float_t dZleft = zPitch - dZright;
+
+ Float_t xpix = nXpix;
+ Float_t dXright = xPitch*(xpix - xPixn);
+ Float_t dXleft = xPitch - dXright;
+
+ Float_t dZprev = 0.;
+ Float_t dZnext = 0.;
+ Float_t dXprev = 0.;
+ Float_t dXnext = 0.;
+
+ for(jz=jzmin; jz <=jzmax; jz++) {
+ if(jz == 1) {
+ dZprev = -zPitch - dZleft;
+ dZnext = -dZleft;
+ }
+ if(jz == 2) {
+ dZprev = -dZleft;
+ dZnext = dZright;
+ }
+ if(jz == 3) {
+ dZprev = dZright;
+ dZnext = dZright + zPitch;
+ }
+ // kz changes from 1 to the fNofPixels(270)
+ Int_t kz = nZpix + jz -2;
+
+ Float_t zArg1 = dZprev/sigmaDif;
+ Float_t zArg2 = dZnext/sigmaDif;
+ Float_t zProb1 = TMath::Erfc(zArg1);
+ Float_t zProb2 = TMath::Erfc(zArg2);
+ Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge;
+
+
+ // ----------- holes diffusion in X(r*phi) direction --------
+
+ if(dZCharge > 1.) {
+ for(jx=jxmin; jx <=jxmax; jx++) {
+ if(jx == 1) {
+ dXprev = -xPitch - dXleft;
+ dXnext = -dXleft;
+ }
+ if(jx == 2) {
+ dXprev = -dXleft;
+ dXnext = dXright;
+ }
+ if(jx == 3) {
+ dXprev = dXright;
+ dXnext = dXright + xPitch;
+ }
+ Int_t kx = nXpix + jx -2;
+
+ Float_t xArg1 = dXprev/sigmaDif;
+ Float_t xArg2 = dXnext/sigmaDif;
+ Float_t xProb1 = TMath::Erfc(xArg1);
+ Float_t xProb2 = TMath::Erfc(xArg2);
+ Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
+
+ if(dXCharge > 1.) {
+ Int_t index = kz-1;
+
+ if (first) {
+ indexRange[0]=indexRange[1]=index;
+ indexRange[2]=indexRange[3]=kx-1;
+ first=kFALSE;
+ }
+
+ indexRange[0]=TMath::Min(indexRange[0],kz-1);
+ indexRange[1]=TMath::Max(indexRange[1],kz-1);
+ indexRange[2]=TMath::Min(indexRange[2],kx-1);
+ indexRange[3]=TMath::Max(indexRange[3],kx-1);
+
+ // build the list of digits for this module
+ Double_t signal=fMapA2->GetSignal(index,kx-1);
+ signal+=dXCharge;
+ fMapA2->SetHit(index,kx-1,(double)signal);
+ } // dXCharge > 1 e-
+ } // jx loop
+ } // dZCharge > 1 e-
+ } // jz loop
+ } // iZi loop
+
+ if (status == 65) { // the step is inside of Si
+ zPix0 = zPix;
+ xPix0 = xPix;
+ }
+ yPrev = yPix;
+
+ if(dray == 0) {
+ GetList(itrack,idhit,pList,indexRange);
+ }
+
+ lasttrack=itrack;
+ } // hit loop inside the module
+
+
+ // introduce the electronics effects and do zero-suppression
+ ChargeToSignal(pList);
+
+ // clean memory
+
+ fMapA2->ClearMap();
+
+
+}
+
+//---------------------------------------------
+void AliITSsimulationSPDdubna::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange)
+{
+ // lop over nonzero digits
+
+
+ //set protection
+ for(int k=0;k<4;k++) {
+ if (indexRange[k] < 0) indexRange[k]=0;
+ }
+
+ for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
+ for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){
+
+ Float_t signal=fMapA2->GetSignal(iz,ix);
+
+ if (!signal) continue;
+
+ Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0!
+ if(!pList[globalIndex]){
+
+ //
+ // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
+ //
+
+ pList[globalIndex] = new Float_t [9];
+
+ // set list to -3
+
+ *pList[globalIndex] = -3.;
+ *(pList[globalIndex]+1) = -3.;
+ *(pList[globalIndex]+2) = -3.;
+ *(pList[globalIndex]+3) = 0.;
+ *(pList[globalIndex]+4) = 0.;
+ *(pList[globalIndex]+5) = 0.;
+ *(pList[globalIndex]+6) = -1.;
+ *(pList[globalIndex]+7) = -1.;
+ *(pList[globalIndex]+8) = -1.;
+
+
+ *pList[globalIndex] = (float)label;
+ *(pList[globalIndex]+3) = signal;
+ *(pList[globalIndex]+6) = (float)idhit;
+ }
+ else{
+
+ // check the signal magnitude
+
+ Float_t highest = *(pList[globalIndex]+3);
+ Float_t middle = *(pList[globalIndex]+4);
+ Float_t lowest = *(pList[globalIndex]+5);
+
+ signal -= (highest+middle+lowest);
+
+ //
+ // compare the new signal with already existing list
+ //
+
+ if(signal<lowest) continue; // neglect this track
+
+ if (signal>highest){
+ *(pList[globalIndex]+5) = middle;
+ *(pList[globalIndex]+4) = highest;
+ *(pList[globalIndex]+3) = signal;
+
+ *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
+ *(pList[globalIndex]+1) = *pList[globalIndex];
+ *pList[globalIndex] = label;
+
+ *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
+ *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
+ *(pList[globalIndex]+6) = idhit;
+ }
+ else if (signal>middle){
+ *(pList[globalIndex]+5) = middle;
+ *(pList[globalIndex]+4) = signal;
+
+ *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
+ *(pList[globalIndex]+1) = label;
+
+ *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
+ *(pList[globalIndex]+7) = idhit;
+ }
+ else{
+ *(pList[globalIndex]+5) = signal;
+ *(pList[globalIndex]+2) = label;
+ *(pList[globalIndex]+8) = idhit;
+ }
+ }
+ } // end of loop pixels in x
+ } // end of loop over pixels in z
+
+
+}
+
+
+//---------------------------------------------
+void AliITSsimulationSPDdubna::ChargeToSignal(Float_t **pList)
+{
+ // add noise and electronics, perform the zero suppression and add the
+ // digit to the list
+
+ AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
+
+
+ Float_t threshold = (float)fResponse->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;iz<fNPixelsZ;iz++){
+ for(Int_t ix=0;ix<fNPixelsX;ix++){
+ electronics = fBaseline + fNoise*gRandom->Gaus();
+ 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;i<fNPixelsZ;i++) {
+ TString spdName("spd_");
+ Char_t pixelz[4];
+ sprintf(pixelz,"%d",i+1);
+ spdName.Append(pixelz);
+ (*fHis)[i] = new TH1F(spdName.Data(),"SPD maps",
+ fNPixelsX,0.,(Float_t) fNPixelsX);
+ }
+}
+
+//____________________________________________
+
+void AliITSsimulationSPDdubna::ResetHistograms()
+{
+ //
+ // Reset histograms for this detector
+ //
+
+ for ( int i=0;i<fNPixelsZ;i++ ) {
+ if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
+ }
+
+}
+
+
+
+
+
+
+
+
+
--- /dev/null
+#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
+
+
+
--- /dev/null
+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: "<<ver<<" ! (must be 5 for the moment)\n";
+ return 12345;
+ }
+
+ if (ver==5) {
+ gROOT->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;
+}
--- /dev/null
+#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 " <<nev<<endl;
+ cout << "nparticles " <<nparticles<<endl;
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+
+ Int_t last_entry=0;
+ timer.Start();
+ ITS->DigitsToRecPoints(nev,last_entry,"All");
+ timer.Stop(); timer.Print();
+ } // event loop
+
+ delete rec0;
+ delete rec1;
+ delete rec2;
+
+ file->Close();
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#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 " <<nev<<endl;
+ if(nev>0) {
+ nparticles = gAlice->GetEvent(nev);
+ gAlice->SetEvent(nev);
+ if(!gAlice->TreeD()) gAlice-> MakeTree("D");
+ ITS->MakeBranch("D");
+ }
+ cout << "nparticles " <<nparticles<<endl;
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+
+ Int_t nbgr_ev=0;
+ if(nsignal) nbgr_ev=Int_t(nev/nsignal);
+ timer.Start();
+ ITS->HitsToDigits(nev,nbgr_ev,size," ","All"," ");
+ timer.Stop(); timer.Print();
+ } // event loop
+
+ delete sim0;
+ delete sim1;
+ delete sim2;
+
+
+ file->Close();
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#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 <<endl;
+ cout << "nparticles " << nparticles <<endl;
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+
+ TTree *TH = gAlice->TreeH();
+ Int_t ntracks = TH->GetEntries();
+ cout<<"ntracks "<<ntracks<<endl;
+
+ Int_t nbytes = 0;
+
+ //AliITSRawClusterSDD *ITSclust;
+ AliITSRawClusterSPD *ITSclust;
+
+// Get pointers to Alice detectors and Digits containers
+ AliITS *ITS = (AliITS*)gAlice->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; mod<nent; mod++) {
+ ITS->ResetClusters();
+ 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;icl<ncl;icl++) {
+ ITSclust = (AliITSRawClusterSPD*)ITSclu->UncheckedAt(icl);
+ printf("%d %d %f %f %f\n",ITSclust->NclZ(),ITSclust->NclX(),ITSclust->Q(),ITSclust->X(),ITSclust->Z());
+ }
+ //for (Int_t icl=0;icl<ncl;icl++) {
+ //ITSclust = (AliITSRawClusterSDD*)ITSclu->UncheckedAt(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 "<<endl;
+
+ file->Close();
+}
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#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 <<endl;
+ cout << "nparticles " << nparticles <<endl;
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+
+ TTree *TH = gAlice->TreeH();
+ Int_t ntracks = TH->GetEntries();
+ cout<<"ntracks "<<ntracks<<endl;
+
+// Get pointers to Alice detectors and Digits containers
+ AliITS *ITS = (AliITS*)gAlice->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; mod<nent; mod++) {
+ AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(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 ="<<mod<<","<<lay<<","<<nclust<<","<<nhits<<endl;
+ for (Int_t clu=0;clu<nclust;clu++) {
+ itsclu = (AliITSRawClusterSPD*)ITSclusters->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;hit<nhits;hit++) {
+
+ // Find coordinate differences between the hit and cluster positions
+ // for the resolution determination.
+
+ itsHit = (AliITShit*)itsModule->GetHit(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 ="<<zhit<<endl;
+ zhit = SPDlength/2 - 10;
+ }
+ if(zhit < 0 && zhit < -SPDlength/2) {
+ //cout<<"!!! z outside ="<<zhit<<endl;
+ zhit = -SPDlength/2 + 10;
+ }
+ if(xhit > SPDwidth/2) {
+ //cout<<"!!! x outside ="<<xhit<<endl;
+ xhit = SPDwidth/2 - 10;
+ }
+ if(xhit < 0 && xhit < -SPDwidth/2) {
+ //cout<<"!!! x outside ="<<xhit<<endl;
+ xhit = -SPDwidth/2 + 10;
+ }
+
+ zhit += SPDlength/2;
+ xhit += SPDwidth/2;
+ Float_t yhit = 10000*itsHit->GetYL();
+
+ if(hitlayer == 1 && hitstat == 66 && yhit > 71) {
+ xhit0 = xhit;
+ zhit0 = zhit;
+ }
+ if(hitlayer == 2 && hitstat == 66 && yhit < -71) {
+ xhit0 = xhit;
+ zhit0 = zhit;
+ }
+
+ if(hitstat != 68) continue; // Take only the hit if the last
+ // track point went out from
+ // the detector.
+
+ if(xhit0 > 9e+4 || zhit0 > 9e+4) continue;
+
+ Float_t xmed = (xhit + xhit0)/2;
+ Float_t zmed = (zhit + zhit0)/2;
+
+ Float_t xdif = xmed - clusterx;
+ Float_t zdif = zmed - clusterz;
+
+
+ // Consider the hits inside of cluster region only
+
+ if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) {
+ icl = 1;
+
+ // part = (TParticle *)particles.UncheckedAt(track);
+ // Int_t partcode = part->GetPdgCode();
+ // Int_t primery = gAlice->GetPrimary(track);
+
+ Int_t parent = itsHit->GetParticle()->GetFirstMother();
+ Int_t partcode = itsHit->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 pmod = itsHit->GetParticle()->P(); // total momentum at the
+ // vertex
+ Float_t energy = itsHit->GetParticle()->Energy(); // energy at the
+ // vertex
+ Float_t mass = itsHit->GetParticle()->GetMass(); // particle mass
+
+ Float_t pz = itsHit->GetParticle()->Pz(); // z momentum componetnt
+ // at the vertex
+ Float_t px = itsHit->GetParticle()->Px(); // z momentum componetnt
+ // at the vertex
+ Float_t py = itsHit->GetParticle()->Py(); // z momentum componetnt
+ // at the vertex
+ Float_t phi = itsHit->GetParticle()->Phi(); // Phi angle at the
+ // vertex
+ Float_t theta = itsHit->GetParticle()->Theta(); // Theta angle at the
+ // vertex
+ //Float_t eta = itsHit->GetParticle()->Eta(); // Pseudo rapidity at the
+ // vertex
+ if((energy-pz) > 0) {
+ Float_t y = 0.5*TMath::Log((energy+pz)/(energy-pz));
+ }else{
+ cout<<" Warning: energy < pz ="<<energy<<","<<pz<<endl;
+ y = 10;
+ }
+ Float_t eta = -TMath::Log(TMath::Tan(theta/2));
+ pmod *= 1.0e+3;
+
+
+ Float_t pxsim = itsHit->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 ="<<occup1<<endl;
+ cout<<" Occupancy for layer-2 ="<<occup2<<endl;
+ // The real occupancy values are:
+ // (for full ALICE event at the full SPD acceptence)
+ // occup1 /= 3932160;
+ // occup2 /= 7864320;
+
+ } // idettype loop
+ } // end if ITS
+} // event loop
+
+
+ // Write and Draw Histogramms and ntuples
+
+
+
+ TFile fhistos("SPD_his_dubna.root","RECREATE");
+
+ ntuple->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"<<endl;
+
+ TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
+ c1->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 "<<endl;
+
+// file->Close();
+}
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff