]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Added in ITS upgrade -
authoramastros <amastros@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Jan 2011 09:47:02 +0000 (09:47 +0000)
committeramastros <amastros@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Jan 2011 09:47:02 +0000 (09:47 +0000)
- constructor for standard ITS configuration
- geometry/segmentation printout summary added.
- example macros to red hits, digits and recpoints

ITS/UPGRADE/AliITSupgrade.cxx
ITS/UPGRADE/AliITSupgrade.h
ITS/UPGRADE/readDigits.C [new file with mode: 0644]
ITS/UPGRADE/readHit.C [new file with mode: 0644]
ITS/UPGRADE/readRecPoint.C [new file with mode: 0644]

index f741b168124cf0ddf364c60f1c1f0b631d74cf76..4e847fab33ce1f21e576b84f7fd18360499daeb9 100644 (file)
@@ -54,6 +54,83 @@ ClassImp(AliITSupgrade)
 {
   //
   // 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):
@@ -64,9 +141,9 @@ AliITSupgrade::AliITSupgrade(const char *name,const char *title, TArrayD widths,
   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),
@@ -84,6 +161,12 @@ AliITSupgrade::AliITSupgrade(const char *name,const char *title, TArrayD widths,
     fRadii.AddAt(radii.At(i),i);
     AliDebug(1,"Creating Volume");
   }
+
+  if(bp){
+    fRadiusBP=radiusBP;
+    fWidthBP=widthBP;
+    fHalfLengthBP=halfLengthBP;
+  }
   Init();
 
 }
@@ -166,6 +249,7 @@ void AliITSupgrade::CreateGeometry()
   TGeoVolumeAssembly *vol= CreateVol();
   gGeoManager->GetVolume("ALIC")->AddNode(vol,0);
   AliInfo("Stop ITS upgrade preliminary version building");
+  PrintSummary();
 } 
 //__________________________________________________________________________________________________
 void AliITSupgrade::Init()
@@ -232,7 +316,7 @@ void AliITSupgrade::StepManager()
     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.
@@ -464,3 +548,64 @@ void AliITSupgrade::SetTreeAddress()
   }
   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;
+ }
+
+}
index 7523bc87e8bebb6b9e27099ee1f5018cd8eb6092..64816d00188b08e4ba1e67ee38a48eda7e78009b 100644 (file)
@@ -36,7 +36,8 @@ class AiITShit;
 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
   
@@ -80,8 +81,15 @@ class AliITSupgrade : public AliITS //TObject-TNamed-AliModule-AliDetector-AliIT
   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;
diff --git a/ITS/UPGRADE/readDigits.C b/ITS/UPGRADE/readDigits.C
new file mode 100644 (file)
index 0000000..4a52c4e
--- /dev/null
@@ -0,0 +1,98 @@
+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();
+  }
+
+
+}
diff --git a/ITS/UPGRADE/readHit.C b/ITS/UPGRADE/readHit.C
new file mode 100644 (file)
index 0000000..9bb5aae
--- /dev/null
@@ -0,0 +1,83 @@
+
+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();
+  }
+
+}
diff --git a/ITS/UPGRADE/readRecPoint.C b/ITS/UPGRADE/readRecPoint.C
new file mode 100644 (file)
index 0000000..9b91dc0
--- /dev/null
@@ -0,0 +1,92 @@
+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();
+  } 
+}
+