#include<TClassTable.h>
#include<TClonesArray.h>
#include<TFile.h>
+#include<TGeoManager.h>
#include<TH1.h>
#include<TH2.h>
+#include <TInterpreter.h>
#include<TObject.h>
#include<TObjArray.h>
#include<TTree.h>
-#include <AliRun.h>
-#include <AliITS.h>
-#include <AliITSgeom.h>
-#include <AliITSDetType.h>
-#include <AliITSRecPoint.h>
-#include <AliITSdigit.h>
-#include <AliITShit.h>
-#include <AliITSmodule.h>
-#include <AliITSsegmentation.h>
-#include <AliITSsegmentationSPD.h>
-#include <AliITSsegmentationSDD.h>
-#include <AliITSsegmentationSSD.h>
-#include <AliRunLoader.h>
-#include <AliITSLoader.h>
-#include <AliHeader.h>
+#include "AliRun.h"
+#include "AliITS.h"
+#include "AliITSgeom.h"
+#include "AliITSDetTypeRec.h"
+#include "AliITSRecPoint.h"
+#include "AliITSRecPoint.h"
+#include "AliITSdigit.h"
+#include "AliITSdigitSSD.h"
+#include "AliITShit.h"
+#include "AliITSmodule.h"
+#include "AliITSsegmentation.h"
+#include "AliITSsegmentationSPD.h"
+#include "AliITSsegmentationSDD.h"
+#include "AliITSsegmentationSSD.h"
+#include "AliRunLoader.h"
+#include "AliITSLoader.h"
+#include "AliHeader.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
#endif
void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb);
Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
+Int_t GetClusCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos);
-Int_t AliITSGeoPlot (Int_t evesel=0, char *opt="All+Rec", char *filename="galice.root", Int_t isfastpoints = 0) {
+Int_t AliITSGeoPlot (Int_t evesel=0, char *opt="All+ClustersV2", TString filename="galice.root", Int_t isfastpoints = 0) {
/*******************************************************************
* This macro displays geometrical information related to the
- * hits, digits and rec points in ITS.
+ * hits, digits and rec points (or V2 clusters) in ITS.
* There are histograms that are not displayed (i.e. energy
* deposition) but that are saved onto a file (see below)
*
* it is wise to redirect the output onto a file
* e.g.: .x AliITSGeoPlot.C("All+Verbose") > out.log
* 3) Rec Points option: "Rec" ---> Uses Rec. Points (default)
- * otherwise ---> uses hits and digits only
+ *
+ * 4) ClustersV2 option: "ClustersV2" ---->Uses ClustersV2
+ * otherwise ---->uses hits and digits only
* Examples:
* .x AliITSGeoPlot(); (All subdetectors; no-verbose; no-recpoints)
* .x AliITSGeoPlot("SPD+SSD+Verbose+Rec");
* isfastpoints: integer. It is set to 0 by defaults. This means that
* slow reconstruction is assumed. If fast recpoint are in
* in use, isfastpoints must be set =1.
- * BEWARE: to analyze digits/recpoints produced with aliroot
- * versions older than v3-07-Release, the input parameters
- * filename, FileDigits, FileRec and isfastpoints must be left
- * to their default value
*
* OUTPUT: It produces a root file with a list of histograms
*
* --- AliITSGeoPlot();
*
* M.Masera 14/05/2001 18:30
- * Last rev. 09/06/2003 17:00 (Adapted to NewIO) m.m.
+ * Last rev. 31/05/2004 14:00 (Clusters V2 added) E.C.
********************************************************************/
//Options
Bool_t All = choice.Contains("All");
Bool_t verbose=choice.Contains("Verbose");
Bool_t userec=choice.Contains("Rec");
+ Bool_t useclustersv2=choice.Contains("ClustersV2");
Int_t retcode=1; //return code
-#if !(!defined(__CINT__) || defined(__MAKECINT__))
if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
+ gInterpreter->ExecuteMacro("loadlibs.C");
}
- else {
-#endif
+ else {
if(gAlice){
- delete gAlice->GetRunLoader();
+ delete AliRunLoader::Instance();
delete gAlice;
gAlice=0;
}
-#if !(!defined(__CINT__) || defined(__MAKECINT__))
}
-#endif
-
-AliRunLoader* rl = AliRunLoader::Open(filename);
- if (rl == 0x0){
- cerr<<"AliITSGeoPlot.C : Can not open session RL=NULL"<< endl;
- return -1;
- }
- Int_t retval = rl->LoadgAlice();
- if (retval){
- cerr<<"AliITSGeoPlot.C : LoadgAlice returned error"<<endl;
- return -1;
- }
- gAlice=rl->GetAliRun();
+ // Set OCDB if needed
+ AliCDBManager* man = AliCDBManager::Instance();
+ if (!man->IsDefaultStorageSet()) {
+ printf("Setting a local default storage and run number 0\n");
+ man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ man->SetRun(0);
+ }
+ else {
+ printf("Using deafult storage \n");
+ }
+ // retrives geometry
+ TString geof(gSystem->DirName(filename));
+ geof += "/geometry.root";
+ TGeoManager::Import(geof.Data());
+ if (!gGeoManager) {
+ cout<<"geometry not found\n";
+ return -1;
+ }
+
+ AliRunLoader* rl = AliRunLoader::Open(filename.Data());
+ if (rl == 0x0){
+ cerr<<"AliITSGeoPlot.C : Can not open session RL=NULL"<< endl;
+ return -1;
+ }
+ Int_t retval = rl->LoadgAlice();
+ if (retval){
+ cerr<<"AliITSGeoPlot.C : LoadgAlice returned error"<<endl;
+ return -1;
+ }
+ gAlice=rl->GetAliRun();
- retval = rl->LoadHeader();
- if (retval){
- cerr<<"AliITSGeoPlot.C : LoadHeader returned error"<<endl;
- return -1;
- }
+ retval = rl->LoadHeader();
+ if (retval){
+ cerr<<"AliITSGeoPlot.C : LoadHeader returned error"<<endl;
+ return -1;
+ }
- AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
+ AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
- if(!ITSloader){
- cerr<<"AliITSGeoPlot.C : ITS loader not found"<<endl;
- return -1;
- }
+ if(!ITSloader){
+ cerr<<"AliITSGeoPlot.C : ITS loader not found"<<endl;
+ return -1;
+ }
- ITSloader->LoadHits("read");
- ITSloader->LoadDigits("read");
- if(isfastpoints==1)ITSloader->SetRecPointsFileName("ITS.FastRecPoints.root");
- ITSloader->LoadRecPoints("read");
- rl->GetEvent(evesel);
- Int_t nparticles = rl->GetHeader()->GetNtrack();
- AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
- ITS->SetTreeAddress();
+ ITSloader->LoadHits("read");
+ ITSloader->LoadDigits("read");
+ if(isfastpoints==1)ITSloader->SetRecPointsFileName("ITS.FastRecPoints.root");
+ ITSloader->LoadRecPoints("read");
+ rl->GetEvent(evesel);
+ Int_t nparticles = rl->GetHeader()->GetNtrack();
+ AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
+ ITS->SetTreeAddress();
if(verbose) {
cout<<" "<<endl<<" "<<endl;
cout<<"******* Event processing started *******"<<endl;
cout<<"Filling modules... It takes a while, now. Please be patient"<<endl;
ITS->FillModules(0,0,nmodules," "," ");
cout<<"ITS modules .... DONE!"<<endl;
+
+ AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
+ detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
+ detTypeRec->SetDefaults();
// DIGITS
TTree *TD = ITSloader->TreeD();
- //RECPOINTS
+ //RECPOINTS (V2 clusters)
TTree *TR = ITSloader->TreeR();
- TClonesArray *ITSrec = ITS->RecPoints();
+ TClonesArray *ITSrec = detTypeRec->RecPoints();
TBranch *branch = 0;
if(userec && TR && ITSrec){
if(isfastpoints==1){
cout<<"WARNING: there are no RECPOINTS on this file ! \n";
cout<<"======================================================= \n \n";
}
+ if(useclustersv2 && TR && ITSrec){
+ branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
+ if(branch)branch->SetAddress(&ITSrec);
+ }
+
+ if(useclustersv2 && (!TR || !ITSrec || !branch)){
+ useclustersv2 = kFALSE;
+ cout<<"\n ======================================================= \n";
+ cout<<"WARNING: there are no CLUSTERSV2 on this file ! \n";
+ cout<<"======================================================= \n \n";
+ }
+
//local variables
Int_t mod; //module number
TH1F *zspd = new TH1F("zspd","Z of digits - SPD",100,-30.,30.);
TH1F *zhspd = new TH1F("zhspd","Z of hits - SPD",100,-30.,30.);
TH1F *zmspd = new TH1F("zmspd","Z of SPD modules",100,-30,30.);
- TH2F *rrspd = new TH2F("rrspd","Radii of recpoints - SPD",50,-10.,10.,50,-10.,10.);
- TH1F *zrspd = new TH1F("zrspd","Z of recpoints - SPD",100,-30.,30.);
+
+ Char_t title1[50]="";
+ Char_t title2[50]="";
+ if(userec){
+ sprintf(title1,"Radii of recpoints - %s","SPD");
+ sprintf(title2,"Z of recpoints - %s","SPD");
+ }
+ if(useclustersv2){
+ sprintf(title1,"Radii of clustersV2 - %s","SPD");
+ sprintf(title2,"Z of clustersV2 - %s","SPD");
+ }
+ TH2F *rrspd = new TH2F("rrspd",title1,50,-10.,10.,50,-10.,10.);
+ TH1F *zrspd = new TH1F("zrspd",title2,100,-30.,30.);
TH1F *enespd = new TH1F("enespd","Energy deposition SPD (KeV)",100,0.,1000.);
histos.AddLast(rspd); // 0
histos.AddLast(rhspd); // 1
TH1F *zsdd = new TH1F("zsdd","Z of digits - SDD",100,-40.,40.);
TH1F *zhsdd = new TH1F("zhsdd","Z of hits - SDD",100,-40.,40.);
TH1F *zmsdd = new TH1F("zmsdd","Z of SDD modules",100,-40,40.);
- TH2F *rrsdd = new TH2F("rrsdd","Radii of recpoints - SDD",50,-40.,40.,50,-40.,40.);
- TH1F *zrsdd = new TH1F("zrsdd","Z of recpoints - SDD",100,-40.,40.);
+ Char_t title3[50];
+ Char_t title4[50];
+ if(userec){
+ sprintf(title3,"Radii of recpoints - %s","SDD");
+ sprintf(title4,"Z of recpoints - %s","SDD");
+ }
+ if(useclustersv2){
+ sprintf(title3,"Radii of clustersV2 - %s","SDD");
+ sprintf(title4,"Z of clustersV2 - %s","SDD");
+ }
+ TH2F *rrsdd = new TH2F("rrsdd",title3,50,-40.,40.,50,-40.,40.);
+ TH1F *zrsdd = new TH1F("zrsdd",title4,100,-40.,40.);
TH1F *enesdd = new TH1F("enesdd","Energy deposition SDD (KeV)",100,0.,1000.);
histos.AddLast(rsdd); // 9
histos.AddLast(rhsdd); // 10
TH1F *zssd = new TH1F("zssd","Z of digits - SSD",100,-70.,70.);
TH1F *zhssd = new TH1F("zhssd","Z of hits - SSD",100,-70.,70.);
TH1F *zmssd = new TH1F("zmssd","Z of SSD modules",100,-70,70.);
- TH2F *rrssd = new TH2F("rrssd","Radii of recpoints - SSD",50,-50.,50.,50,-50.,50.);
- TH1F *zrssd = new TH1F("zrssd","Z of recpoints - SSD",100,-70.,70.);
+ Char_t title5[50];
+ Char_t title6[50];
+ if(userec){
+ sprintf(title5,"Radii of recpoints - %s","SSD");
+ sprintf(title6,"Z of recpoints - %s","SSD");
+ }
+ if(useclustersv2){
+ sprintf(title5,"Radii of clustersV2 - %s","SSD");
+ sprintf(title6,"Z of clustersV2 - %s","SSD");
+ }
+
+ TH2F *rrssd = new TH2F("rrssd",title5,50,-50.,50.,50,-50.,50.);
+ TH1F *zrssd = new TH1F("zrssd",title6,100,-70.,70.);
TH1F *enessd = new TH1F("enessd","Energy deposition SSD (KeV)",100,0.,1000.);
histos.AddLast(rssd); // 18
histos.AddLast(rhssd); // 19
//
// Loop on subdetectors
//
+ cout<<"CALL GETITSGEOM \n \n \n";
AliITSgeom *geom = ITS->GetITSgeom();
TString detna; // subdetector name
for(Int_t subd=0;subd<3;subd++){
cout<<"======================================================= \n \n";
}
// Get segmentation model
- AliITSDetType *iDetType=ITS->DetType(subd);
if(subd==0)detna="SPD";
if(subd==1)detna="SDD";
if(subd==2)detna="SSD";
- AliITSsegmentation *seg=(AliITSsegmentation*)iDetType->GetSegmentationModel();
+ AliITSsegmentation *seg=(AliITSsegmentation*)detTypeRec->GetSegmentationModel(subd);
// Loop on modules
first = geom->GetStartDet(subd);
last = geom->GetLastDet(subd);
//RecPoints
if(userec){
- ITS->ResetRecPoints();
- // TR->GetEvent(mod);
+ detTypeRec->ResetRecPoints();
branch->GetEvent(mod);
TH2F *bidi=(TH2F*)histos.At(6+subd*9);
TH1F *uni=(TH1F*)histos.At(7+subd*9);
nrecp=GetRecCoor(geom,ITSrec,mod,bidi,uni,verbose);
}
+ if(useclustersv2){
+ detTypeRec->ResetRecPoints();
+ branch->GetEvent(mod);
+ TH2F *bidi=(TH2F*)histos.At(6+subd*9);
+ TH1F *uni=(TH1F*)histos.At(7+subd*9);
+ nrecp=GetClusCoor(geom,ITSrec,mod,bidi,uni,verbose);
+
+ }
// Digits
if(usedigits){
- ITS->ResetDigits();
+ detTypeRec->ResetDigits();
nbytes += TD->GetEvent(mod);
GetDigits(seg,geom,ITSdigits,subd,mod,verbose,histos);
}
h1tmp=(TH1F*)histos.At(4+subd*9);
h1tmp->Draw();
- if(userec){
+ if(userec || useclustersv2){
current->cd(5);
h2tmp=(TH2F*)histos.At(6+9*subd);
h2tmp->Draw();
}
+Int_t GetClusCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
+
+ AliITSgeom *geom = (AliITSgeom*)ge;
+ Int_t nrecp = ITSrec->GetEntries();
+ if(nrecp>0){
+ Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
+ Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
+ if(verb){
+ cout<<"-------------------------------------------------------"<<endl;
+ cout<<"Number of CLUSTERS for module "<<mod<<": "<<nrecp<<endl;
+ }
+ for(Int_t irec=0;irec<nrecp;irec++) {
+ AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
+ Double_t rot[9];
+ geom->GetRotMatrix(mod,rot);
+ Int_t lay,lad,det;
+ geom->GetModuleId(mod,lay,lad,det);
+ Float_t tx,ty,tz;
+ geom->GetTrans(lay,lad,det,tx,ty,tz);
+
+ Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
+ Double_t phi1=TMath::Pi()/2+alpha;
+ if(lay==1) phi1+=TMath::Pi();
+
+ Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
+ Float_t r=tx*cp+ty*sp;
+ gc[0]= r*cp - recp->GetY()*sp;
+ gc[1]= r*sp + recp->GetY()*cp;
+ gc[2]=recp->GetZ();
+
+ if(verb){
+ Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
+ cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" ";
+ cout<<r<<endl;
+ cout<<"Associated track "<<recp->GetLabel(0)<<endl;
+ }
+ h2->Fill(gc[0],gc[1]);
+ h1->Fill(gc[2]);
+
+ }
+ }
+ return nrecp;
+}
Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
AliITSgeom *geom = (AliITSgeom*)ge;
Int_t nrecp = ITSrec->GetEntries();
}
for(Int_t irec=0;irec<nrecp;irec++) {
AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
- lc[0]=recp->GetX();
- lc[2]=recp->GetZ();
+ lc[0]=recp->GetDetLocalX();
+ lc[2]=recp->GetDetLocalZ();
geom->LtoG(mod,lc,gc);
if(verb){
cout<<"recp # "<<irec<<" local coordinates. lx= "<<lc[0]<<" lz= ";
}
for (Int_t digit=0;digit<ndigits;digit++) {
digs = (AliITSdigit*)ITSdigits->UncheckedAt(digit);
- Int_t iz=digs->fCoord1; // cell number z
- Int_t ix=digs->fCoord2; // cell number x
+ Int_t iz=digs->GetCoord1(); // cell number z
+ Int_t ix=digs->GetCoord2(); // cell number x
// Get local coordinates of the element
if(subd<2){
seg->DetToLocal(ix,iz,lcoor[0],lcoor[2]);
for(Int_t digi2=0;digi2<ndigits;digi2++){
if(ssdone[digi2]==0 && impaired){
AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(digi2);
- if(dig2->fCoord1 != iz && dig2->GetTrack(0)==digs->GetTrack(0) && dig2->GetTrack(0)>0){
+ if(dig2->GetCoord1() != iz && dig2->GetTrack(0)==digs->GetTrack(0) && dig2->GetTrack(0)>0){
ssdone[digi2]=2;
pair[digit]=digi2;
- if(pside)nstrip=dig2->fCoord2;
- if(!pside)pstrip=dig2->fCoord2;
+ if(pside)nstrip=dig2->GetCoord2();
+ if(!pside)pstrip=dig2->GetCoord2();
impaired=kFALSE;
}
}
}
}
if(subd<2 || (subd==2 && ssdone[digit]==1)){
- Int_t coor1=digs->fCoord1;
- Int_t coor2=digs->fCoord2;
+ Int_t coor1=digs->GetCoord1();
+ Int_t coor2=digs->GetCoord2();
Int_t tra0=digs->GetTrack(0);
if(verbose){
cout<<"digit # "<<digit<<" fCoord1= "<<coor1<<" fCoord2= ";
cout<<endl;
Int_t dtmp=pair[digit];
AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(dtmp);
- Int_t coor1b=dig2->fCoord1;
- Int_t coor2b=dig2->fCoord2;
+ Int_t coor1b=dig2->GetCoord1();
+ Int_t coor2b=dig2->GetCoord2();
Int_t tra0b=dig2->GetTrack(0);
cout<<"(digit paired with digit #"<<dtmp<<endl;
cout<<"with fCoord1= "<<coor1b<<" fCoord2= "<<coor2b;
geom->LtoG(mod,lcoor,gcoor); // global coord. in cm
ragdig=sqrt(gcoor[0]*gcoor[0]+gcoor[1]*gcoor[1]);
if(verbose){
- cout<<"global coordinates "<<gcoor[0]<<" "<<gcoor[1];
- cout<<" "<<gcoor[2]<<" Radius "<<ragdig<<endl;
+ cout<<"global coordinates "<<gcoor[0]<<" "<<gcoor[1];
+ cout<<" "<<gcoor[2]<<" Radius "<<ragdig<<endl;
}
//Fill histograms
TH2F *bidi = (TH2F*)histos.At(subd*9);
} // loop on digits for this module
} // if(ndigits>0....
}
+