{
//
// default ctor
+ //
+}
+
+//__________________________________________________________________________________________________
+AliITSupgrade::AliITSupgrade(const char *name,const char *title, Bool_t isBeamPipe):
+ AliITS(name,title),
+ fWidths(0),
+ fRadii(0),
+ fRadiiCu(0),
+ fWidthsCu(0),
+ fCopper(0),
+ fBeampipe(isBeamPipe),
+ fRadiusBP(0),
+ fWidthBP(0),
+ fHalfLengthBP(0),
+ fNlayers(0),
+ fHalfLength(0),
+ fSdigits(0),
+ fDigits(0),
+ fSegmentation(0x0)
+{
+ fNlayers = 6;
+ fWidths.Set(fNlayers);
+ fRadii.Set(fNlayers);
+ fRadiiCu.Set(fNlayers);
+ fWidthsCu.Set(fNlayers);
+ fCopper.Set(fNlayers);
+ fHalfLength.Set(fNlayers);
+
+ // Default values are used in order to simulate the standard ITS material budget (see The ALICE Collaboration et al 2008 JINST 3 S08002 - Figure 3.2).
+ // The cell segmentation is chosen to achieve the tracking resolution as described in Table 3.2 of The ALICE Collaboration et al 2008 JINST 3 S08002,
+ // apart from SPD0 where the 9 um value is considered : resolution*sqrt(12)
+
+ Double_t xsz[6]={31.18*1e-04,41.6*1e-04,131.6*1e-04,131.6*1e-04,69.3*1e-04,69.3*1e-04};
+ Double_t zsz[6]={416*1e-04,416*1e-04,97*1e-04,97*1e-04,2875*1e-04,2875*1e-04};
+ Double_t Lhalf[6]={80.,80.,80.,80.,80.,80.};
+ Double_t r[6]={4.,7.6,14.9,23.8,39.1,43.6};
+ Double_t thick[6]={75.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04};
+
+ Int_t nlayers =6;
+ TArrayD xsizeSi(nlayers);
+ TArrayD zsizeSi(nlayers);
+
+ // adding beam pipe upgrade
+ if(isBeamPipe){
+ fRadiusBP = 2.9;
+ fWidthBP = 0.08;
+ fHalfLengthBP = Lhalf[0];
+ }
+
+ // setting the geonetry parameters and the segmentation
+ Int_t npixHalf[6], npixR[6];
+ Double_t HalfL[6], xszInt[6];
+ Int_t c[6]={1,1,1,1,1,1};
+
+ for(Int_t i=0; i<nlayers; i++){
+ // recalc length
+ npixHalf[i]=(Int_t)(Lhalf[i]/zsz[i]);
+ HalfL[i]=npixHalf[i]*zsz[i];
+ // recalc segmentation
+ npixR[i] = (Int_t)(2*TMath::Pi()*r[i]/xsz[i]);
+ xszInt[i]= 2*TMath::Pi()*r[i]/npixR[i];
+ xsizeSi.AddAt(xszInt[i],i);
+ zsizeSi.AddAt(zsz[i],i);
+
+ fHalfLength.AddAt(HalfL[i],i);
+
+ fRadii.AddAt(r[i],i);
+ fWidths.AddAt(thick[i],i);
+
+ fRadiiCu.AddAt(r[i]+thick[i],i);
+ fWidthsCu.AddAt(0.015,i);//micron
+ fCopper.AddAt(c[i],i);
+ }
+
+ SetFullSegmentation(xsizeSi,zsizeSi);
+ Init();
}
//__________________________________________________________________________________________________
AliITSupgrade::AliITSupgrade(const char *name,const char *title, TArrayD widths, TArrayD radii,TArrayD halfLengths, TArrayD radiiCu, TArrayD widthsCu, TArrayS copper,Bool_t bp,Double_t radiusBP, Double_t widthBP, Double_t halfLengthBP):
fWidthsCu(widthsCu),
fCopper(copper),
fBeampipe(bp),
- fRadiusBP(radiusBP),
- fWidthBP(widthBP),
- fHalfLengthBP(halfLengthBP),
+ fRadiusBP(0),
+ fWidthBP(0),
+ fHalfLengthBP(0),
fNlayers(widths.GetSize()),
fHalfLength(halfLengths),
fSdigits(0),
fRadii.AddAt(radii.At(i),i);
AliDebug(1,"Creating Volume");
}
+
+ if(bp){
+ fRadiusBP=radiusBP;
+ fWidthBP=widthBP;
+ fHalfLengthBP=halfLengthBP;
+ }
Init();
}
TGeoVolumeAssembly *vol= CreateVol();
gGeoManager->GetVolume("ALIC")->AddNode(vol,0);
AliInfo("Stop ITS upgrade preliminary version building");
+ PrintSummary();
}
//__________________________________________________________________________________________________
void AliITSupgrade::Init()
hit.SetStartStatus(status);
return; // don't save entering hit.
}
- // Fill hit structure with this new hit.
+ // Fill hit structure with this new hit.
new((*fHits)[fNhits++]) AliITShit(hit); // Use Copy Construtor.
// Save old position... for next hit.
}
AliDebug(1,"Stop.");
}
+//____________________________________________________________________________________________________
+void AliITSupgrade::PrintSummary()
+{
+
+ // pdg booklet or http://pdg.lbl.gov/2010/reviews/rpp2010-rev-atomic-nuclear-prop.pdf
+
+ // X0 as X0[g/cm^2]/density
+ Double_t beX0= 65.19/1.848;
+ Double_t cuX0= 12.86/8.96;
+ Double_t siX0= 21.82/2.329;
+
+ if(!fSegmentation) fSegmentation= new AliITSsegmentationUpgrade();
+ printf(" %10s %10s %10s %10s %10s %10s %10s %10s\n","Name","R [cm]","x/X0 [%]","Phi res[um]","Z res[um]","X segm[um]","Y segm[um]","widths [um]");
+
+ printf(" Beampipe %10.3f %10f %10s %10s %10s %10s %10.1f\n", fRadiusBP, 100.*fWidthBP/beX0,"-","-","-","-",fWidthBP*1.e+4);
+ //else printf(" ------- No Beam Pipe Upgrade ---------------- \n");
+ for(Int_t i=0; i< fNlayers; i++){
+ printf(" Si Layer %1i %10.3f %10f %10.2f %10.2f %10.2f %10.2f %10.1f\n",i, fRadii.At(i), 100.*fWidths.At(i)/siX0, (fSegmentation->GetCellSizeX(i) *1.e+4)/TMath::Sqrt(12.), (fSegmentation->GetCellSizeZ(i)*1.e+4)/TMath::Sqrt(12.),fSegmentation->GetCellSizeX(i) *1.e+4, fSegmentation->GetCellSizeZ(i)*1.e+4,fWidths.At(i)*1.e+4);
+ printf(" Cu Layer %10.3f %10f %10s %10s %10s %10s %10.1f\n", fRadiiCu.At(i), 100.*fWidthsCu.At(i)/cuX0,"-","-","-","-",fWidthsCu.At(i)*1.e+4);
+ }
+}
+//____________________________________________________________________________________________________
+void AliITSupgrade::SetSegmentationX(Double_t x, Int_t lay){
+TFile *f = TFile::Open("Segmentation.root","UPDATE");
+if(!f) {
+ AliError("\n\n\n Segmentation.root file does not exist. The segmentation is 0! \n\n\n");
+ return;
+ }
+else {
+ TArrayD *a = (TArrayD*)f->Get("CellSizeX");
+ if(!a) {
+ AliError("CessSizeX array does not exist!");
+ return;
+ }
+ a->AddAt(x,lay);
+ f->WriteObjectAny(a,"TArrayD","CellSizeX");
+ f->Close();
+ delete f;
+ }
+}
+//____________________________________________________________________________________________________
+void AliITSupgrade::SetSegmentationZ(Double_t z, Int_t lay){
+
+TFile *f = TFile::Open("Segmentation.root","UPDATE");
+if(!f) {
+ AliError("\n\n\n Segmentation.root file does not exist. The segmentation is 0! \n\n\n");
+ return;
+ }
+else {
+ TArrayD *a = (TArrayD*)f->Get("CellSizeZ");
+ if(!a) {
+ AliError("CellSizeZ array does not exist!");
+ return;
+ }
+ a->AddAt(z,lay);
+ f->WriteObjectAny(a,"TArrayD","CellSizeZ");
+ f->Close();
+ delete f;
+ }
+
+}
class AliITSupgrade : public AliITS //TObject-TNamed-AliModule-AliDetector-AliITS-AliITSupgrade
{
public:
- AliITSupgrade(); //default ctor
+ AliITSupgrade(); //default ctor
+ AliITSupgrade(const char *name, const char *title, Bool_t isBeamPipe=kTRUE); //ctor for standard ITS
AliITSupgrade(const char *name, const char *title, TArrayD widths, TArrayD radii,TArrayD halfLengths, TArrayD radiiCu, TArrayD widthsCu, TArrayS copper,Bool_t bp, Double_t radiusBP, Double_t widthPB, Double_t halfLengthsBP); //ctor
virtual ~AliITSupgrade(); //dtor
TGeoVolumeAssembly * CreateVol();
void SetFullSegmentation(TArrayD xsize, TArrayD zsize);
+ void SetRadius(Double_t r, Int_t lay) {fRadii.AddAt(r,lay);}
+ void SetWidth(Double_t w, Int_t lay) {fWidths.AddAt(w,lay);}
+ void SetRadiusCu(Double_t rCu, Int_t lay) {fRadiiCu.AddAt(rCu,lay);}
+ void SetWidthCu(Double_t wCu, Int_t lay) {fWidthsCu.AddAt(wCu,lay);}
+ void SetSegmentationX(Double_t x, Int_t lay);
+ void SetSegmentationZ(Double_t z, Int_t lay);
void StepHistory();
+ void PrintSummary();
protected:
TArrayD fWidths;
--- /dev/null
+void readDigits(){
+
+ gSystem->Load("libITSUpgradeBase");
+ gSystem->Load("libITSUpgradeSim");
+ gROOT->SetStyle("Plain");
+
+ const Int_t nLayers = 6;
+ Int_t nbins=100;
+ Int_t xmin=0;
+ Int_t xmax=50000;//00*1e-09;
+
+
+ TH1D *hNelSDig[nLayers], *hNelDig[nLayers];
+ for(Int_t i=0; i< nLayers; i++ ) {
+ hNelSDig[i] = new TH1D(Form("hNelSDig%i",i),Form("electron distribution in SDigits [ Layer %i] ",i),nbins,xmin,xmax);
+ hNelSDig[i]->SetXTitle("N electrons");
+ hNelDig[i] = new TH1D(Form("hNelDig%i",i),Form("electron distribution in Digits [ Layer %i] ",i),nbins,xmin,xmax);
+ hNelDig[i]->SetXTitle("N electrons");
+ }
+
+ const Int_t layers = 6;
+
+ gAlice=NULL;
+ AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
+ runLoader->LoadgAlice();
+
+ gAlice = runLoader->GetAliRun();
+
+ runLoader->LoadHeader();
+ runLoader->LoadKinematics();
+ runLoader->LoadSDigits();
+ runLoader->LoadDigits();
+
+ AliLoader *dl = runLoader->GetDetectorLoader("ITS");
+
+
+ //SDIGITS INIT
+ TTree * sDigTree = 0x0;
+ TObjArray *sDigList= new TObjArray(layers);
+ for(Int_t i=0; i<layers; i++) sDigList->AddAt(new TClonesArray("AliITSDigitUpgrade"),i);
+
+ //DIGITS INIT
+ TTree * digTree = 0x0;
+ TObjArray *digList= new TObjArray(layers);
+ for(Int_t i=0; i<layers; i++) digList->AddAt(new TClonesArray("AliITSDigitUpgrade"),i);
+
+ printf("N Events : %i \n",(Int_t)runLoader->GetNumberOfEvents());
+
+ for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
+ printf("\n Event %i \n",iEvent);
+ runLoader->GetEvent(iEvent);
+ AliStack *s = runLoader->Stack();
+ sDigTree=dl->TreeS();
+ digTree=dl->TreeD();
+
+ for(Int_t ilayer=0; ilayer<layers; ilayer ++) {
+ digTree->SetBranchAddress(Form("Layer%d",ilayer),&(*digList)[ilayer]);
+ sDigTree->SetBranchAddress(Form("Layer%d",ilayer),&(*sDigList)[ilayer]);
+ }
+
+ // LOOPS
+
+ digTree->GetEntry(0);
+ sDigTree->GetEntry(0);
+
+ for(Int_t il=0; il<layers; il++){
+ TClonesArray *pArrDig= (TClonesArray*)digList->At(il);
+ TClonesArray *pArrSdig= (TClonesArray*)sDigList->At(il);
+ for(Int_t isdentry=0; isdentry< pArrSdig->GetEntries(); isdentry++) {
+ AliITSDigitUpgrade *pSdig = (AliITSDigitUpgrade*)pArrSdig->At(isdentry);
+ cout << " --- SDigit generated by a "<<s->Particle(pSdig->GetTrack(0))->GetName()<<"----"<<endl;
+ hNelSDig[pSdig->GetLayer()]->Fill(pSdig->GetNelectrons());
+ pSdig->PrintInfo();
+ }
+
+ for(Int_t ientry=0; ientry< pArrDig->GetEntries(); ientry++) {
+ cout << " --- Digit ---"<<endl;
+ AliITSDigitUpgrade *pDig = (AliITSDigitUpgrade*)pArrDig->At(ientry);
+ hNelDig[pDig->GetLayer()]->Fill(pDig->GetNelectrons());
+ pDig->PrintInfo();
+ }
+ }
+
+ }//event loop
+ TCanvas *cSd = new TCanvas("cSd","Summable Digits",1000,800);
+ cSd->Divide(3,2);
+ TCanvas *cD = new TCanvas("cD","Digits",1000,800);
+ cD->Divide(3,2);
+
+ for(Int_t ip =1; ip<=6; ip++){
+ cSd->cd(ip);
+ hNelSDig[ip-1]->Draw();
+ cD->cd(ip);
+ hNelDig[ip-1]->Draw();
+ }
+
+
+}
--- /dev/null
+
+void readHit(){
+ gROOT->SetStyle("Plain");
+
+ gSystem->Load("libITSUpgradeBase");
+ gSystem->Load("libITSUpgradeSim");
+
+ TH2F *xyGlob = new TH2F("xyGlob"," X - Y Global coordinates ",100,-50,50,100,-50,50);
+ xyGlob->SetXTitle("cm");
+ xyGlob->SetMarkerStyle(7);
+ TH1F *zGlob = new TH1F("zGlob", " Z Global coordinates ",200, -100,100 );
+ zGlob->SetXTitle("cm");
+
+ Int_t nbins=100;
+ Int_t xmin=0;
+ Int_t xmax=0.1;//00*1e-09;
+
+ const Int_t nLayers = 6;
+
+ TH1D *hDeLoss[nLayers];
+ for(Int_t i=0; i< nLayers; i++ ) {
+ hDeLoss[i] = new TH1D(Form("hDeLossl%i",i),Form("E loss distribution [ Layer %i] ",i),nbins,xmin,xmax);
+ hDeLoss[i]->SetXTitle("GeV");
+ }
+
+ gAlice=NULL;
+ AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
+ runLoader->LoadgAlice();
+
+ gAlice = runLoader->GetAliRun();
+
+ runLoader->LoadHeader();
+ runLoader->LoadKinematics();
+ runLoader->LoadSDigits();
+ runLoader->LoadHits();
+
+ AliLoader *dl = runLoader->GetDetectorLoader("ITS");
+
+ //HIT INIT
+
+ TTree *hitTree = 0x0;
+ TClonesArray *hitList=new TClonesArray("AliITShit");
+
+ AliITSsegmentationUpgrade *segmentation = new AliITSsegmentationUpgrade();
+ for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
+
+ runLoader->GetEvent(iEvent);
+
+ hitTree=dl->TreeH();
+ hitTree->SetBranchAddress("ITSupgrade",&hitList);
+ for(Int_t iEnt=0;iEnt<hitTree->GetEntries();iEnt++){//entries loop degli hits
+ hitTree->GetEntry(iEnt);
+ for(Int_t iHit=0; iHit<hitList->GetEntries();iHit++){
+
+ AliITShit *pHit = (AliITShit*)hitList->At(iHit);
+
+ if(pHit->GetParticle()->IsPrimary()){
+ Double_t xg,yg,zg=0.;
+ pHit->GetPositionG(xg,yg,zg);
+ xyGlob->Fill(xg,yg);
+ zGlob->Fill(zg);
+ hDeLoss[pHit->GetModule()]->Fill(pHit->GetIonization());
+ } // is primary
+ }//loop hit
+ }//entryloopHitList
+
+ }//event loop
+
+ TCanvas *xyCanv = new TCanvas("xyCanv","Hit X-Y positions",500,500);
+ xyCanv->cd();
+ xyGlob->Draw();
+ TCanvas *zCanv = new TCanvas("zCanv","Hit Z positions",500,500);
+ zCanv->cd();
+ zGlob->Draw();
+
+ TCanvas *c = new TCanvas("c","E loss distribution per layer",1000,800);
+ c->Divide(3,2);
+ for(Int_t ip =1; ip<=6; ip++){
+ c->cd(ip);
+ hDeLoss[ip-1]->Draw();
+ }
+
+}
--- /dev/null
+void readRecPoint(){
+ gSystem->Load("libITSUpgradeSim");
+ gSystem->Load("libITSUpgradeBase");
+ gSystem->Load("libITSUpgradeRec");
+ gROOT->SetStyle("Plain");
+ Int_t nbins=100;
+ Int_t xmin=0;
+ Int_t xmax=50000;//00*1e-09;
+
+ const Int_t nLayers = 6;
+
+ AliITSsegmentationUpgrade *seg = new AliITSsegmentationUpgrade();
+ if(!seg){
+ printf("no segmentation info available... Exiting");
+ return;
+ }
+
+ TH1D *hNel[nLayers];
+ for(Int_t i=0; i< nLayers; i++ ) {
+ hNel[i] = new TH1D(Form("hNel%i",i),Form("cluster charge distribution [ Layer %i] ",i),nbins,xmin,xmax);
+ hNel[i]->SetXTitle("N electrons");
+ }
+ TH1D * type = new TH1D("hCluType"," cluster type" , 50,0,15 );
+
+ TH2F *xyGlob = new TH2F("xyGlob"," X - Y Global coordinates ",100,-50,50,100,-50,50);
+ xyGlob->SetXTitle("cm");
+ xyGlob->SetMarkerStyle(7);
+ TH1F *zGlob = new TH1F("zGlob", " Z Global coordinates ",200, -100,100 );
+ zGlob->SetXTitle("cm");
+
+
+ gAlice=NULL;
+ AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
+ runLoader->LoadgAlice();
+
+ gAlice = runLoader->GetAliRun();
+
+ runLoader->LoadHeader();
+ runLoader->LoadKinematics();
+ runLoader->LoadRecPoints();
+
+ AliITSLoader *dl = (AliITSLoader*)runLoader->GetDetectorLoader("ITS");
+
+
+ TTree *clusTree = 0x0;
+
+ TClonesArray statITSCluster("AliITSRecPoint");
+
+ for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
+ runLoader->GetEvent(iEvent);
+ clusTree=dl->TreeR();
+ TClonesArray *ITSCluster = &statITSCluster;
+ TBranch* itsClusterBranch=clusTree->GetBranch("ITSRecPoints");
+ if (!itsClusterBranch) {
+ printf("can't get the branch with the ITS clusters ! \n");
+ return;
+ }
+ itsClusterBranch->SetAddress(&ITSCluster);
+ clusTree->GetEntry(0);
+ Double_t charge=0.;
+ Int_t nCluster = ITSCluster->GetEntriesFast();
+ for(Int_t i=0; i<nCluster; i++){
+ AliITSRecPoint *recp = (AliITSRecPoint*)ITSCluster->UncheckedAt(i);
+ Double_t xyz[3]={-1,-1,-1};
+ seg->DetToGlobal(recp->GetLayer(), recp->GetDetLocalX(), recp->GetDetLocalZ(), xyz[0],xyz[1],xyz[2]) ;
+ xyGlob->Fill(xyz[0],xyz[1]);
+ zGlob->Fill(xyz[2]);
+ charge=recp->GetQ();
+ // cout<< "layer "<< recp->GetLayer() << " local system X "<< recp->GetDetLocalX() << " Z "<< recp->GetDetLocalZ() <<endl;
+ type->Fill(recp->GetType());
+ hNel[recp->GetLayer()]->Fill(charge);
+ }
+
+ }
+
+ TCanvas *xyCanv = new TCanvas("xvCanvClus","RecPoint X-Y Positions",500,500);
+ xyCanv->cd();
+ xyGlob->Draw();
+ TCanvas *zCanv = new TCanvas("zCanvClus","RecPoint Z Positions",500,500);
+ zCanv->cd();
+ zGlob->Draw();
+ new TCanvas();
+ type->Draw();
+
+ TCanvas *c = new TCanvas("c","Cluster charge distribution",1000,800);
+ c->Divide(3,2);
+ for(Int_t ip =1; ip<=6; ip++){
+ c->cd(ip);
+ hNel[ip-1]->Draw();
+ }
+}
+