/*
$Log$
+Revision 1.25 2000/11/12 22:38:05 barbera
+Added header file for the SPD Bari model
+
+Revision 1.24 2000/10/09 22:18:12 barbera
+Bug fixes from MAriana to le AliITStest.C run correctly
+
+Revision 1.23 2000/10/05 20:47:42 nilsen
+fixed dependencies of include files. Tryed but failed to get a root automaticly
+generates streamer function to work. Modified SetDefaults.
+
+Revision 1.9.2.15 2000/10/04 16:56:40 nilsen
+Needed to include stdlib.h
+
+=======
+Revision 1.22 2000/10/04 19:45:52 barbera
+Corrected by F. Carminati for v3.04
+
Revision 1.21 2000/10/02 21:28:08 fca
Removal of useless dependecies via forward declarations
// futher information.
//
///////////////////////////////////////////////////////////////////////////////
-
+#include <stdlib.h>
#include <TMath.h>
#include <TRandom.h>
#include <TBranch.h>
#include "AliITSDetType.h"
#include "AliITSClusterFinder.h"
#include "AliITSsimulation.h"
+#include "AliITSresponse.h"
#include "AliITSsegmentationSPD.h"
#include "AliITSresponseSPD.h"
+#include "AliITSresponseSPDbari.h"
#include "AliITSsegmentationSDD.h"
#include "AliITSresponseSDD.h"
#include "AliITSsegmentationSSD.h"
#include "AliMC.h"
#include "stdlib.h"
+#include "AliITStrack.h"
+#include "AliITSiotrack.h"
+#include "AliITStracking.h"
+#include "../TPC/AliTPC.h"
+#include "../TPC/AliTPCParam.h"
+
const Int_t AliITS::fgkNTYPES=3;
ClassImp(AliITS)
delete fHits;
delete fDigits;
delete fRecPoints;
- if(fIdName!=0) delete[] fIdName;
+// delete fIdName; // TObjArray of TObjStrings
+ if(fIdName!=0) delete[] fIdName; // Array of TStrings
if(fIdSens!=0) delete[] fIdSens;
if(fITSmodules!=0) {
this->ClearModules();
// and sets the default segmentation, response, digit and raw cluster classes
// Therefore it should be called after a call to CreateGeometry.
//
-
-
- SetDefaults();
-
Int_t i;
+
cout << endl;
for(i=0;i<30;i++) cout << "*";cout << " ITS_INIT ";
for(i=0;i<30;i++) cout << "*";cout << endl;
+//
+ SetDefaults();
+// TObjArray of TObjStrings
+// for(i=0;i<fIdN;i++) fIdSens[i] = gMC->VolId((fIdName->At(i))->GetName());
+// Array of TStrings
for(i=0;i<fIdN;i++) fIdSens[i] = gMC->VolId(fIdName[i]);
- //
+//
for(i=0;i<70;i++) cout << "*";
cout << endl;
}
{
// sets the default segmentation, response, digit and raw cluster classes
- cout << "AliITS::SetDefaults" << endl;
+ printf("SetDefaults\n");
AliITSDetType *iDetType;
+
//SPD
- AliITSsegmentationSPD *seg0=new AliITSsegmentationSPD(fITSgeom);
- AliITSresponseSPD *resp0=new AliITSresponseSPD();
iDetType=DetType(0);
- if (!iDetType->GetSegmentationModel()) SetSegmentationModel(0,seg0);
- if (!iDetType->GetResponseModel()) SetResponseModel(0,resp0);
+ if (!iDetType->GetSegmentationModel()) {
+ AliITSsegmentationSPD *seg0=new AliITSsegmentationSPD(fITSgeom);
+ SetSegmentationModel(0,seg0);
+ }
+ if (!iDetType->GetResponseModel()) {
+ SetResponseModel(0,new AliITSresponseSPD());
+ }
// set digit and raw cluster classes to be used
- const char *kData0=resp0->DataType();
+
+ const char *kData0=(iDetType->GetResponseModel())->DataType();
if (strstr(kData0,"real")) {
iDetType->ClassNames("AliITSdigit","AliITSRawClusterSPD");
} else iDetType->ClassNames("AliITSdigitSPD","AliITSRawClusterSPD");
// SDD //
- AliITSresponseSDD *resp1=new AliITSresponseSDD();
- AliITSsegmentationSDD *seg1=new AliITSsegmentationSDD(fITSgeom,resp1);
iDetType=DetType(1);
- //printf("SetDefaults: iDetType %p\n",iDetType);
- if (!iDetType->GetSegmentationModel()) SetSegmentationModel(1,seg1);
- //printf("SetDefaults: segm %p\n",iDetType->GetSegmentationModel());
- if (!iDetType->GetResponseModel()) SetResponseModel(1,resp1);
- //printf("SetDefaults: resp %p\n",iDetType->GetResponseModel());
- const char *kData1=resp1->DataType();
- const char *kopt=resp1->ZeroSuppOption();
+ if (!iDetType->GetResponseModel()) {
+ SetResponseModel(1,new AliITSresponseSDD());
+ }
+ AliITSresponse *resp1=iDetType->GetResponseModel();
+ if (!iDetType->GetSegmentationModel()) {
+ AliITSsegmentationSDD *seg1=new AliITSsegmentationSDD(fITSgeom,resp1);
+ SetSegmentationModel(1,seg1);
+ }
+ const char *kData1=(iDetType->GetResponseModel())->DataType();
+ const char *kopt=iDetType->GetResponseModel()->ZeroSuppOption();
if ((!strstr(kopt,"2D")) && (!strstr(kopt,"1D")) || strstr(kData1,"real") ) {
iDetType->ClassNames("AliITSdigit","AliITSRawClusterSDD");
} else iDetType->ClassNames("AliITSdigitSDD","AliITSRawClusterSDD");
// SSD
- AliITSsegmentationSSD *seg2=new AliITSsegmentationSSD(fITSgeom);
- AliITSresponseSSD *resp2=new AliITSresponseSSD();
iDetType=DetType(2);
- if (!iDetType->GetSegmentationModel()) SetSegmentationModel(2,seg2);
- if (!iDetType->GetResponseModel()) SetResponseModel(2,resp2);
- const char *kData2=resp2->DataType();
+ if (!iDetType->GetSegmentationModel()) {
+ AliITSsegmentationSSD *seg2=new AliITSsegmentationSSD(fITSgeom);
+ SetSegmentationModel(2,seg2);
+ }
+ if (!iDetType->GetResponseModel()) {
+ SetResponseModel(2,new AliITSresponseSSD());
+ }
+ const char *kData2=(iDetType->GetResponseModel())->DataType();
if (strstr(kData2,"real")) {
iDetType->ClassNames("AliITSdigit","AliITSRawClusterSSD");
} else iDetType->ClassNames("AliITSdigitSSD","AliITSRawClusterSSD");
}
}
+
+
+
//_____________________________________________________________________________
void AliITS::SetDefaultSimulation()
{
AliDetector::Streamer(R__b);
R__b >> fIdN;
R__b.ReadArray(fIdSens);
+ if(fIdName!=0) delete[] fIdName; // Array of TStrings
+ fIdName = new TString[fIdN];
for(i=0;i<fIdN;i++) fIdName[i].Streamer(R__b);
R__b >> fITSgeom;
R__b >> fITSmodules;
}
+//________________________________________________________________
+
+AliITStrack AliITS::Tracking(AliITStrack &track, AliITStrack *reference,TObjArray *fastpoints, Int_t
+**vettid, Bool_t flagvert ) {
+
+//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
+
+
+ TList *list= new TList();
+
+ AliITStrack tr(track);
+
+ list->AddLast(&tr);
+
+ Double_t Pt=(tr).GetPt();
+ //cout << "\n Pt = " << Pt <<"\n";
+
+
+ AliITStracking obj(list, reference, this, fastpoints,TMath::Abs(Pt),vettid, flagvert);
+ list->Delete();
+ delete list;
+
+ Int_t itot=-1;
+ TVector VecTotLabref(18);
+ for(Int_t lay=5; lay>=0; lay--) {
+ TVector VecLabref(3);
+ VecLabref=(*reference).GetLabTrack(lay);
+ for(Int_t k=0; k<3; k++) {itot++; VecTotLabref(itot)=VecLabref(k);}
+ }
+ Long_t labref;
+ Int_t freq;
+ (*reference).Search(VecTotLabref, labref, freq);
+
+ if(freq < 5) labref=-labref;
+ (*reference).SetLabel(labref);
+
+ return *reference;
+
+}
+
+
+
+//________________________________________________________________
+
+
+void AliITS::DoTracking(Int_t evNumber, Int_t min_t, Int_t max_t, TFile *file, Bool_t flagvert) {
+
+// ex macro for tracking ITS
+
+ printf("begin DoTracking - file %p\n",file);
+
+ const char *pname="75x40_100x60";
+
+ struct GoodTrack {
+ Int_t lab,code;
+ Float_t px,py,pz,x,y,z,pxg,pyg,pzg,ptg;
+ Bool_t flag;
+ };
+
+
+ gAlice->GetEvent(0);
+
+ AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
+ AliTPCParam *digp = (AliTPCParam*)file->Get(pname);
+ if (digp!=0) TPC->SetParam(digp);
+
+ GoodTrack gt[7000];
+ Int_t ngood=0;
+ ifstream in("good_tracks");
+
+ cerr<<"Reading good tracks...\n";
+ while (in>>gt[ngood].lab>>gt[ngood].code
+ >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
+ >>gt[ngood].x >>gt[ngood].y >>gt[ngood].z
+ >>gt[ngood].pxg >>gt[ngood].pyg >>gt[ngood].pzg
+ >>gt[ngood].ptg >>gt[ngood].flag) {
+ ngood++;
+ cerr<<ngood<<'\r';
+ if (ngood==7000) {
+ cerr<<"Too many good tracks !\n";
+ break;
+ }
+ }
+ if (!in.eof()) cerr<<"Read error (good_tracks) !\n";
+
+
+// Load tracks
+ TFile *tf=TFile::Open("tpctracks.root");
+ if (!tf->IsOpen()) {cerr<<"Can't open tpctracks.root !\n"; return ;}
+ TObjArray tracks(200000);
+ TTree *tracktree=(TTree*)tf->Get("TreeT");
+ TBranch *tbranch=tracktree->GetBranch("tracks");
+ Int_t nentr=(Int_t)tracktree->GetEntries();
+ for (Int_t k=0; k<nentr; k++) {
+ AliTPCtrack *iotrack=new AliTPCtrack;
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(k);
+ tracks.AddLast(iotrack);
+ }
+ tf->Close();
+
+
+ Int_t nt = tracks.GetEntriesFast();
+ cerr<<"Number of found tracks "<<nt<<endl;
+
+ TVector DataOut(9);
+ Int_t kkk=0;
+
+ Double_t ptg=0.,pxg=0.,pyg=0.,pzg=0.;
+
+ ////////////////////////////// good tracks definition in TPC ////////////////////////////////
+
+ ofstream out1 ("AliITSTrag.out");
+ for (Int_t i=0; i<ngood; i++) out1 << gt[i].ptg << "\n";
+ out1.close();
+
+
+ TVector vec(5);
+ TTree *TR=gAlice->TreeR();
+ Int_t nent=(Int_t)TR->GetEntries();
+ TClonesArray *recPoints = RecPoints();
+ Int_t numbpoints;
+ Int_t totalpoints=0;
+ Int_t *np = new Int_t[nent];
+ Int_t **vettid = new Int_t* [nent];
+ for (Int_t mod=0; mod<nent; mod++) {
+ vettid[mod]=0;
+ this->ResetRecPoints();
+ gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
+ numbpoints = recPoints->GetEntries();
+ totalpoints+=numbpoints;
+ np[mod] = numbpoints;
+ //cout<<" mod = "<<mod<<" numbpoints = "<<numbpoints<<"\n"; getchar();
+ vettid[mod] = new Int_t[numbpoints];
+ for (Int_t i=0;i<numbpoints; i++) *(vettid[mod]+i)=0;
+ }
+
+ AliTPCtrack *track;
+
+
+ if(min_t < 0) {min_t = 0; max_t = nt-1;}
+
+/*
+ ///////////////////////////////// Definition of vertex end its error ////////////////////////////
+ ////////////////////////// In the future it will be given by a method ///////////////////////////
+ Double_t Vx=0.;
+ Double_t Vy=0.;
+ Double_t Vz=0.;
+
+ Float_t sigmavx=0.0050; // 50 microns
+ Float_t sigmavy=0.0050; // 50 microns
+ Float_t sigmavz=0.010; // 100 microns
+
+ //Vx+=gRandom->Gaus(0,sigmavx); Vy+=gRandom->Gaus(0,sigmavy); Vz+=gRandom->Gaus(0,sigmavz);
+ TVector vertex(3), ervertex(3)
+ vertex(0)=Vx; vertex(1)=Vy; vertex(2)=Vz;
+ ervertex(0)=sigmavx; ervertex(1)=sigmavy; ervertex(2)=sigmavz;
+ /////////////////////////////////////////////////////////////////////////////////////////////////
+*/
+
+ //TDirectory *savedir=gDirectory;
+
+ TTree tracktree1("TreeT","Tree with ITS tracks");
+ AliITSiotrack *iotrack=0;
+ tracktree1.Branch("ITStracks","AliITSiotrack",&iotrack,32000,0);
+
+ ofstream out ("AliITSTra.out");
+
+ for (int j=min_t; j<=max_t; j++) {
+ track=(AliTPCtrack*)tracks.UncheckedAt(j);
+ Int_t flaglab=0;
+ if (!track) continue;
+ ////// elimination of not good tracks ////////////
+ Int_t ilab=TMath::Abs(track->GetLabel());
+
+ for (Int_t i=0;i<ngood;i++) {
+ //cout<<" ilab, gt[i].lab = "<<ilab<<" "<<gt[i].lab<<"\n"; getchar();
+ if (ilab==gt[i].lab) {
+ flaglab=1;
+ ptg=gt[i].ptg;
+ pxg=gt[i].pxg;
+ pyg=gt[i].pyg;
+ pzg=gt[i].pzg;
+ break;
+ }
+ }
+ //cout<<" j flaglab = " <<j<<" "<<flaglab<<"\n"; getchar();
+ if (!flaglab) continue;
+ //cout<<" j = " <<j<<"\n"; getchar();
+ /*
+ ////// old propagation to the end of TPC //////////////
+ Double_t xk=76.;
+ track->PropagateTo(xk);
+ xk-=0.11;
+ track->PropagateTo(xk,42.7,2.27); //C
+ xk-=2.6;
+ track->PropagateTo(xk,36.2,1.98e-3); //C02
+ xk-=0.051;
+ track->PropagateTo(xk,42.7,2.27); //C
+ ///////////////////////////////////////////////////
+ */
+
+ ////// new propagation to the end of TPC //////////////
+ Double_t xk=77.415;
+ track->PropagateTo(xk, 28.94, 1.204e-3); //Ne
+ xk -=0.01;
+ track->PropagateTo(xk, 44.77, 1.71); //Tedlar
+ xk -=0.04;
+ track->PropagateTo(xk, 44.86, 1.45); //Kevlar
+ xk -=2.0;
+ track->PropagateTo(xk, 41.28, 0.029); //Nomex
+ xk-=16;
+ track->PropagateTo(xk,36.2,1.98e-3); //C02
+ xk -=0.01;
+ track->PropagateTo(xk, 24.01, 2.7); //Al
+ xk -=0.01;
+ track->PropagateTo(xk, 44.77, 1.71); //Tedlar
+ xk -=0.04;
+ track->PropagateTo(xk, 44.86, 1.45); //Kevlar
+ xk -=0.5;
+ track->PropagateTo(xk, 41.28, 0.029); //Nomex
+
+ ///////////////////////////////////////////////////////////////
+
+ ///////////////////////////////////////////////////////////////
+ AliITStrack trackITS(*track);
+ AliITStrack result(*track);
+ AliITStrack primarytrack(*track);
+
+///////////////////////////////////////////////////////////////////////////////////////////////
+ TVector Vgeant(3);
+ Vgeant=result.GetVertex();
+
+ // Definition of Dv and Zv for vertex constraint
+ // Double_t sigmaDv=0.0050; Double_t sigmaZv=0.010;
+ Double_t sigmaDv=0.0015; Double_t sigmaZv=0.0015;
+ Double_t uniform= gRandom->Uniform();
+ Double_t signdv;
+ if(uniform<=0.5) signdv=-1.;
+ else
+ signdv=1.;
+
+ Double_t Vr=TMath::Sqrt(Vgeant(0)*Vgeant(0)+ Vgeant(1)*Vgeant(1));
+ Double_t Dv=gRandom->Gaus(signdv*Vr,(Float_t)sigmaDv);
+ Double_t Zv=gRandom->Gaus(Vgeant(2),(Float_t)sigmaZv);
+
+ //cout<<" Dv e Zv = "<<Dv<<" "<<Zv<<"\n";
+ trackITS.SetDv(Dv); trackITS.SetZv(Zv);
+ trackITS.SetsigmaDv(sigmaDv); trackITS.SetsigmaZv(sigmaZv);
+ result.SetDv(Dv); result.SetZv(Zv);
+ result.SetsigmaDv(sigmaDv); result.SetsigmaZv(sigmaZv);
+ primarytrack.SetDv(Dv); primarytrack.SetZv(Zv);
+ primarytrack.SetsigmaDv(sigmaDv); primarytrack.SetsigmaZv(sigmaZv);
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+
+ primarytrack.PrimaryTrack();
+ TVector d2=primarytrack.Getd2();
+ TVector tgl2=primarytrack.Gettgl2();
+ TVector dtgl=primarytrack.Getdtgl();
+ trackITS.Setd2(d2); trackITS.Settgl2(tgl2); trackITS.Setdtgl(dtgl);
+ result.Setd2(d2); result.Settgl2(tgl2); result.Setdtgl(dtgl);
+ /*
+ trackITS.SetVertex(vertex); trackITS.SetErrorVertex(ervertex);
+ result.SetVertex(vertex); result.SetErrorVertex(ervertex);
+ */
+ Tracking(trackITS,&result,recPoints,vettid, flagvert);
+
+ cout<<" progressive track number = "<<j<<"\r";
+ // cout<<" progressive track number = "<<j<<"\n";
+ Long_t labITS=result.GetLabel();
+ // cout << " ITS track label = " << labITS << "\n";
+ int lab=track->GetLabel();
+ //cout << " TPC track label = " << lab <<"\n";
+ //result.GetClusters(); getchar(); //to print the cluster matrix
+
+//propagation to vertex
+
+ Double_t rbeam=3.;
+ result.Propagation(rbeam);
+ TMatrix *cov;
+ cov=&result.GetCMatrix();
+ Double_t pt=TMath::Abs(result.GetPt());
+ Double_t Dr=result.GetD();
+ Double_t Z=result.GetZ();
+ Double_t tgl=result.GetTgl();
+ Double_t C=result.GetC();
+ Double_t Cy=C/2.;
+ Double_t Dz=Z-(tgl/Cy)*TMath::ASin(result.arga(rbeam));
+ Dz-=Vgeant(2);
+ // cout<<" Dr e dz alla fine = "<<Dr<<" "<<Dz<<"\n"; getchar();
+ Double_t phi=result.Getphi();
+ Double_t phivertex = phi - TMath::ASin(result.argA(rbeam));
+ Double_t duepi=2.*TMath::Pi();
+ if(phivertex>duepi) phivertex-=duepi;
+ if(phivertex<0.) phivertex+=duepi;
+ Double_t Dtot=TMath::Sqrt(Dr*Dr+Dz*Dz);
+//////////////////////////////////////////////////////////////////////////////////////////
+ Int_t NumofCluster, idmodule,idpoint;
+ NumofCluster=result.GetNumClust();
+ if(NumofCluster >=5) {
+
+
+ AliITSiotrack outtrack;
+
+ iotrack=&outtrack;
+
+ iotrack->SetStatePhi(phi);
+ iotrack->SetStateZ(Z);
+ iotrack->SetStateD(Dr);
+ iotrack->SetStateTgl(tgl);
+ iotrack->SetStateC(C);
+ Double_t radius=result.Getrtrack();
+ iotrack->SetRadius(radius);
+ Int_t charge;
+ if(C>0.) charge=-1; else charge=1;
+ iotrack->SetCharge(charge);
+
+
+
+ iotrack->SetCovMatrix(cov);
+
+ Double_t px=pt*TMath::Cos(phi);
+ Double_t py=pt*TMath::Sin(phi);
+ Double_t pz=pt*tgl;
+
+ Double_t xtrack=Dr*TMath::Sin(phi);
+ Double_t ytrack=Dr*TMath::Cos(phi);
+ Double_t ztrack=Dz+Vgeant(2);
+
+
+ iotrack->SetPx(px);
+ iotrack->SetPy(py);
+ iotrack->SetPz(pz);
+ iotrack->SetX(xtrack);
+ iotrack->SetY(ytrack);
+ iotrack->SetZ(ztrack);
+ iotrack->SetLabel(labITS);
+
+ for( Int_t il=0;il<6; il++){
+ iotrack->SetIdPoint(il,result.GetIdPoint(il));
+ iotrack->SetIdModule(il,result.GetIdModule(il));
+ }
+ tracktree1.Fill();
+
+ //cout<<" labITS = "<<labITS<<"\n";
+ //cout<<" phi z Dr tgl C = "<<phi<<" "<<Z<<" "<<Dr<<" "<<tgl<<" "<<C<<"\n"; getchar();
+
+ DataOut(kkk) = ptg; kkk++; DataOut(kkk)=labITS; kkk++; DataOut(kkk)=lab; kkk++;
+
+ for (Int_t il=0;il<6;il++) {
+ idpoint=result.GetIdPoint(il);
+ idmodule=result.GetIdModule(il);
+ *(vettid[idmodule]+idpoint)=1;
+ iotrack->SetIdPoint(il,idpoint);
+ iotrack->SetIdModule(il,idmodule);
+ }
+
+ Double_t difpt= (pt-ptg)/ptg*100.;
+ DataOut(kkk)=difpt; kkk++;
+ Double_t lambdag=TMath::ATan(pzg/ptg);
+ Double_t lam=TMath::ATan(tgl);
+ Double_t diflam = (lam - lambdag)*1000.;
+ DataOut(kkk) = diflam; kkk++;
+ Double_t phig=TMath::ATan(pyg/pxg);
+ Double_t phi=phivertex;
+ Double_t difphi = (phi - phig)*1000.;
+ DataOut(kkk)=difphi; kkk++;
+ DataOut(kkk)=Dtot*1.e4; kkk++;
+ DataOut(kkk)=Dr*1.e4; kkk++;
+ DataOut(kkk)=Dz*1.e4; kkk++;
+ for (int r=0; r<9; r++) { out<<DataOut(r)<<" ";}
+ out<<"\n";
+ kkk=0;
+
+
+ } // end if on NumofCluster
+ //gObjectTable->Print(); // stampa memoria
+ } // end for (int j=min_t; j<=max_t; j++)
+
+ out.close();
+
+
+ static Bool_t first=kTRUE;
+ static TFile *tfile;
+
+ if(first) {
+ tfile=new TFile("itstracks.root","RECREATE");
+ //cout<<"I have opened itstracks.root file "<<endl;
+ }
+ first=kFALSE;
+ tfile->cd();
+ tfile->ls();
+
+ char hname[30];
+ sprintf(hname,"TreeT%d",evNumber);
+
+ tracktree1.Write(hname);
+
+
+
+ TTree *fAli=gAlice->TreeK();
+ TFile *fileAli=0;
+
+ if (fAli) fileAli =fAli->GetCurrentFile();
+ fileAli->cd();
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////
+
+ printf("delete vectors\n");
+ if(np) delete [] np;
+ if(vettid) delete [] vettid;
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff