--- /dev/null
+rm galice.root
+rm AliTPCclusters.root
+rm AliTPCtracks.root
+rm AliTPCtracksSorted.root
+rm AliITSclustersV2.root
+rm AliITStracksV2.root
+rm AliITStracksV2Pid.root
+rm itstracks.root
+rm good_tracks_tpc
+rm itsgood_tracks
--- /dev/null
+TFile* AccessFile(TString inFile="galice.root", TString acctype="R");
+void writeAR(TFile * fin, TFile *fou);
+
+Int_t AliITSHits2SDigitsDubna(Int_t evNumber1=0,Int_t evNumber2=0,
+ TString inFile ="galice.root",
+ TString outFile="galiceS.root"){
+
+ // Dynamically link some shared libs
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ } // end if
+
+ // Connect the Root Galice file containing Geometry, Kine and Hits
+
+ TFile *file;
+ if(outFile.Data() == inFile.Data()){
+ file = AccessFile(inFile,"U");
+ }
+ else {
+ file = AccessFile(inFile);
+ }
+
+ TFile *file2 = 0; // possible output file for TreeS
+
+ if(!(outFile.Data() == inFile.Data())){
+ // open output file and create TreeS on it
+ file2 = gAlice->InitTreeFile("S",outFile);
+ }
+
+ AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
+ if (!ITS) {
+ cerr<<"AliITSHits2DigitsDefault.C : AliITS object not found on file"
+ << endl;
+ return 3;
+ } // end if !ITS
+ if(!(ITS->GetITSgeom())){
+ cerr << " AliITSgeom not found. Can't digitize with out it." << endl;
+ return 4;
+ } // end if
+
+ // For old files, must change SPD noise.
+ AliITSresponseSPDdubna *resp0 = new AliITSresponseSPDdubna();
+ if(ITS->DetType(0)->GetResponseModel() !=0){
+ delete ((AliITSresponse*)ITS->DetType(0)->GetResponseModel());
+ ITS->DetType(0)->ResponseModel(0);
+ } // end if
+ ITS->DetType(0)->ResponseModel(resp0);
+ AliITSsegmentationSPD *seg0 = (AliITSsegmentationSPD*)ITS->DetType(0)->GetSegmentationModel();
+ AliITSsimulationSPDdubna *sim0 = new AliITSsimulationSPDdubna(seg0,resp0);
+ if(ITS->DetType(0)->GetSimulationModel() !=0){
+ delete ((AliITSsimulation*)ITS->DetType(0)->GetSimulationModel());
+ ITS->DetType(0)->SimulationModel(0);
+ } // end if
+ ITS->DetType(0)->SimulationModel(sim0);
+
+ TStopwatch timer;
+ timer.Start();
+ for(Int_t nevent = evNumber1; nevent <= evNumber2; nevent++){
+ gAlice->GetEvent(nevent);
+ if(!gAlice->TreeS() && file2 == 0){
+ cout << "Having to create the SDigits Tree." << endl;
+ gAlice->MakeTree("S");
+ } // end if !gAlice->TreeS()
+ if(file2)gAlice->MakeTree("S",file2);
+ // make branch
+ ITS->MakeBranch("S");
+ ITS->SetTreeAddress();
+ cout<<"Making ITS SDigits for event "<<nevent<<endl;
+ TStopwatch timer;
+ Long_t size0 = file->GetSize();
+ ITS->Hits2SDigits();
+ }
+ timer.Stop();
+ timer.Print();
+
+ // write the AliRun object to the output file
+ if(file2)writeAR(file,file2);
+
+ delete gAlice; gAlice=0;
+ file->Close();
+}
+
+//-------------------------------------------------------------------
+TFile * AccessFile(TString FileName, TString acctype){
+
+ // Function used to open the input file and fetch the AliRun object
+
+ TFile *retfil = 0;
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(FileName);
+ if (file) {file->Close(); delete file; file = 0;}
+ if(acctype.Contains("U")){
+ file = new TFile(FileName,"update");
+ }
+ if(acctype.Contains("N") && !file){
+ file = new TFile(FileName,"recreate");
+ }
+ if(!file) file = new TFile(FileName); // default readonly
+ if (!file->IsOpen()) {
+ cerr<<"Can't open "<<FileName<<" !" << endl;
+ return retfil;
+ }
+
+ // Get AliRun object from file or return if not on file
+ if (gAlice) {delete gAlice; gAlice = 0;}
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (!gAlice) {
+ cerr << "AliRun object not found on file"<< endl;
+ return retfil;
+ }
+ return file;
+}
+
+//-------------------------------------------------------------------
+void writeAR(TFile * fin, TFile *fou) {
+ TDirectory *current = gDirectory;
+ TTree *Te;
+ TTree *TeNew;
+ AliHeader *alhe = new AliHeader();
+ Te = (TTree*)fin->Get("TE");
+ Te->SetBranchAddress("Header",&alhe);
+ Te->SetBranchStatus("*",1);
+ fou->cd();
+ TeNew = Te->CloneTree();
+ TeNew->Write(0,TObject::kOverwrite);
+ gAlice->Write(0,TObject::kOverwrite);
+ current->cd();
+ delete alhe;
+ cout<<"AliRun object written to file\n";
+}
+
+
+
+
+
--- /dev/null
+/****************************************************************************
+ * *
+ * This macro determines the impact parameter of each track in pp event *
+ * with respect to the primaty vertex position reconstructed using all *
+ * the other tracks in the event. This procedure is necessary to avoid *
+ * biases in the impact parameter estimate. *
+ * *
+ * Output: TNtuple (ntd0) on file ImpactParameters.root *
+ * There is one entry per track. The elements of the ntuple are: *
+ * 1) Event: Event number *
+ * 2) nTraks: Number of tracks for this event *
+ * 3) pt: Transverse momentum for the current track *
+ * 4) d0rphi: transverse impact parameter of the current track *
+ * 5) vx: x coordinate of the vertex *
+ * 6) vy: y coordinate of the vertex *
+ * origin: A.Dainese andrea.dainese@pd.infn.it *
+ ****************************************************************************/
+#if !defined(__CINT__) || defined(__MAKECINT__)
+//-- --- standard headers-------------
+#include <Riostream.h>
+//--------Root headers ---------------
+#include <TSystem.h>
+#include <TFile.h>
+#include <TNtuple.h>
+#include <TString.h>
+#include <TStopwatch.h>
+#include <TObject.h>
+#include <TVector3.h>
+#include <TTree.h>
+#include <TParticle.h>
+#include <TArray.h>
+//----- AliRoot headers ---------------
+#include "alles.h"
+#include "AliRun.h"
+#include "AliKalmanTrack.h"
+#include "AliITStrackV2.h"
+#include "AliHeader.h"
+#include "AliGenEventHeader.h"
+#include "AliV0vertex.h"
+#include "AliV0vertexer.h"
+#include "AliITSVertex.h"
+#include "AliITSVertexer.h"
+#include "AliITSVertexerTracks.h"
+#endif
+//-------------------------------------
+
+// field (T)
+const Double_t kBz = 0.4;
+
+// this function ckecks file existence
+Bool_t GetInputFile();
+void MoveTracks(TTree& itsTree,TObjArray& array);
+
+
+void AliITSImpactParametersPP(Int_t evFirst=0,Int_t evLast=0) {
+
+ const Char_t *name="AliITSImpactParametersPP";
+ cerr<<'\n'<<name<<"...\n";
+ gBenchmark->Start(name);
+
+ AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
+
+
+ // check existence of input file
+ if(!GetInputFile()) { cerr<<"No tracks file found"<<endl; return; }
+
+ TString outName("ImpactParameters.root");
+
+ // primary vertex
+ Double_t v1[3] = {0.,0.,0,};
+
+ Int_t ev,itsEntries,i;
+ Int_t nTotEv=0;
+ Double_t pt,d0rphi;
+ Char_t trksName[100];
+ AliITStrackV2 *itstrack = 0;
+
+ // output ntuple
+ TNtuple *ntd0 = new TNtuple("ntd0","ntd0","Event:nTrks:pt:d0rphi:vx:vy");
+
+ // create the AliITSVertexerTracks object
+ AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
+ vertexer1->SetMinTracks(3);
+ vertexer1->SetDebug(0);
+ Int_t skipped[2];
+
+ // Open file with ITS tracks
+ TFile* itstrks = TFile::Open("AliITStracksV2.root");
+
+
+ // loop on events in file
+ for(ev=evFirst; ev<=evLast; ev++) {
+ printf(" --- Processing event %d ---\n",ev);
+ sprintf(trksName,"TreeT_ITS_%d",ev);
+
+ // tracks from ITS
+ TTree *itsTree =(TTree*)itstrks->Get(trksName);
+ if(!itsTree) continue;
+ itsEntries = (Int_t)itsTree->GetEntries();
+ printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
+
+ TObjArray tArray(itsEntries);
+ MoveTracks(*itsTree,tArray);
+
+ // count the total number of events
+ nTotEv++;
+
+ // loop on tracks
+ for(i=0; i<itsEntries; i++) {
+ itstrack = (AliITStrackV2*)tArray.At(i);
+
+ // pt
+ pt = 1./TMath::Abs(itstrack->Get1Pt());
+
+ // primary vertex from other tracks in event
+ Bool_t goodVtx = kFALSE;
+ vertexer1->SetVtxStart(0.,0.);
+ skipped[0] = i;
+ skipped[1] = -1;
+ vertexer1->SetSkipTracks(1,skipped);
+ AliITSVertex *vertex1 =
+ (AliITSVertex*)vertexer1->VertexOnTheFly(*itsTree);
+ if(vertex1->GetNContributors()>0) goodVtx = kTRUE;
+ vertex1->GetXYZ(v1);
+ //vertex1->PrintStatus();
+ // impact parameter w.r.t. v1 (in microns)
+ d0rphi = 10000.*itstrack->GetD(v1[0],v1[1]);
+ delete vertex1;
+
+ // fill ntuple
+ if (goodVtx) ntd0->Fill(ev,itsEntries,pt,d0rphi,v1[0],v1[1]);
+
+ itstrack = 0;
+ } // loop on tracks
+
+ delete itsTree;
+
+ } // loop on events in file
+
+
+ printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
+
+ delete vertexer1;
+
+ itstrks->Close();
+
+ // store ntuple in a file
+ TFile* outroot = new TFile(outName.Data(),"recreate");
+ ntd0->Write();
+ outroot->Close();
+ delete outroot;
+
+
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+
+ return;
+}
+//___________________________________________________________________________
+Bool_t GetInputFile() {
+ TString itsName("AliITStracksV2.root");
+
+ if(gSystem->AccessPathName(itsName.Data(),kFileExists)) return kFALSE;
+
+ return kTRUE;
+}
+//___________________________________________________________________________
+void MoveTracks(TTree& itsTree,TObjArray& array) {
+//
+// this function creates two TObjArrays with tracks
+//
+
+ Int_t entr = (Int_t)itsTree.GetEntries();
+
+ // trasfer tracks from tree to array
+ for(Int_t i=0; i<entr; i++) {
+
+ AliITStrackV2 *t = new AliITStrackV2;
+ itsTree.SetBranchAddress("tracks",&t);
+
+ itsTree.GetEvent(i);
+
+ array.AddLast(t);
+
+ }
+
+ return;
+}
+
+
+
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "iostream.h"
+#include "TDatetime.h"
+#include "STEER/AliRun.h"
+#include "STEER/AliRunDigitizer.h"
+#include "ITS/AliITSDigitizer.h"
+#include "ITS/AliITS.h"
+#include "ITS/AliITSDetType.h"
+#include "ITS/AliITSresponseSDD.h"
+#include "TStopwatch.h"
+
+#endif
+
+TFile* AccessFile(TString inFile="galice.root", TString acctype="R");
+void writeAR(TFile * fin, TFile *fou);
+Int_t ChangeITSDefaults(TFile *hitfile,AliITS *ITS,TString opt="");
+//#define DEBUG
+Int_t AliITSMerge(TString digFile="galiceMD.root",
+ TString sdigFileSig="galiceS_sig.root",
+ TString sdigFileBg="",TString opt=""){
+ // Standeard ITS SDigits to Digits, with posible merging.
+
+ // Dynamically link some shared libs
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ } // end if
+ // Connect the Root Galice file containing Geometry, Kine and Hits
+ if (gAlice) {delete gAlice; gAlice = 0;}
+
+ TFile *sdigfilesig = 0; // pointer to signal input file.
+ TFile *sdigfilebg = 0; // pointer to background input file.
+ TFile *digfile = 0; // possible output file for TreeD
+
+ // Setup to copy gAlice and the event tree to the output file.
+
+ if(digFile.CompareTo(sdigFileSig)==0){//write output to same file as input.
+ sdigfilesig = AccessFile(sdigFileSig,"U");//input file open for update.
+ }else{ // different output file then input file.
+ sdigfilesig = AccessFile(sdigFileSig,"R");//input file open read only
+ digfile = new TFile(digFile,"NEW");
+ } // end if digFile == hitFile.
+
+ AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
+ if (!ITS) {
+ cerr<<"AliITSHits2DigitsDefault.C : AliITS object not found on file"
+ << endl;
+ return 3;
+ } // end if !ITS
+
+ ChangeITSDefaults(sdigfilesig,ITS,opt);
+ // write the AliRun object to the output file if different from input file.
+ if(digfile){
+ writeAR(sdigfilesig,digfile);
+ digfile->Close(); // Manager will open in update mode.
+ digfile = 0;
+ } // end if digfile
+ sdigfilesig->Close();
+ sdigfilesig = 0;
+// delete gAlice; gAlice=0; // there is a problem with deleting gAlice????
+
+ AliRunDigitizer *manager;
+ if(sdigFileBg.CompareTo("")==0) { // do not merge.
+ manager = new AliRunDigitizer(1,1);
+ }else{
+ manager = new AliRunDigitizer(2,1);
+ manager->SetInputStream(0,sdigFileSig.Data());
+ manager->SetInputStream(1,sdigFileBg.Data());
+ } // end if
+ if (digFile.CompareTo(sdigFileSig) !=0) {
+ manager->SetOutputFile(digFile);
+ } // end if
+// manager->SetCopyTreesFromInput(0);
+ AliITSDigitizer *dITS = new AliITSDigitizer(manager);
+ if(opt.Contains("ROI")){
+ cout << "Region of Interest selected" << endl;
+ dITS->SetByRegionOfInterestFlag(1);
+ }else{
+ cout << "Digizing everthing" << endl;
+ dITS->SetByRegionOfInterestFlag(0);
+ } // end if
+
+ TStopwatch timer;
+ timer.Start();
+ manager->Exec("all");
+ timer.Stop();
+ timer.Print();
+
+ if(digfile!=0){
+ cout << digFile << " size =" << digfile->GetSize() << endl;
+ }else if(sdigfilesig!=0){
+ cout << sdigFileSig << " size =" << sdigfilesig->GetSize() << endl;
+ } // end if sdigfile!=0
+
+
+ if(digfile) digfile->Close();
+ if(sdigfilesig) sdigfilesig->Close();
+ if(sdigfilebg) sdigfilebg->Close();
+
+ // There is a problem deleting gAlice. It is related to the destructor
+ // and the fTreeE in AliRun. So for now it is not deleted.
+// delete gAlice; // digfile is closed by deleting gAlice if != hitfile.
+// gAlice = 0;
+ delete manager;
+}
+//______________________________________________________________________
+TFile * AccessFile(TString FileName, TString acctype){
+ // Function used to open the input file and fetch the AliRun object
+
+ TFile *retfil = 0;
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(FileName);
+ if(file) {
+ file->Close();
+ delete file;
+ file = 0;
+ } // end if file
+ if(acctype.Contains("U")){
+ file = new TFile(FileName,"UPDATE");
+ } // end if open for update
+ if(acctype.Contains("N") && !file){
+ file = new TFile(FileName,"RECREATE");
+ } // end if open a new file
+ if(!file) file = new TFile(FileName,"READ"); // default readonly
+ if (!file->IsOpen()) {
+ cerr << "Can't open " << FileName << " !" << endl;
+ return retfil;
+ } // end if error opeing file
+
+ // Get AliRun object from file or return if not on file
+ if (gAlice) {delete gAlice; gAlice = 0;}
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (!gAlice) {
+ cerr << "AliRun object not found on file "<< FileName << "!" << endl;
+ file->Close(); // close file and return error.
+ return retfil;
+ } // end if !gAlice
+ return file;
+}
+//______________________________________________________________________
+void writeAR(TFile * fin, TFile *fou) {
+ TDirectory *current = gDirectory;
+ TTree *TeOld;
+ TTree *TeNew;
+ AliHeader *alhe = new AliHeader();
+ TeOld = (TTree*)fin->Get("TE");
+ TeOld->SetBranchAddress("Header",&alhe);
+ TeOld->SetBranchStatus("*",1);
+ fou->cd();
+ TeNew = TeOld->CloneTree();
+ TeNew->Write(0,TObject::kOverwrite);
+ gAlice->Write(0,TObject::kOverwrite);
+ current->cd();
+ delete alhe;
+#ifdef DEBUG
+ cout << "AliRun object written to file" << endl;
+#endif
+}
+//______________________________________________________________________
+Int_t ChangeITSDefaults(TFile *hitfile,AliITS *ITS,TString opt){
+
+ TDatime *ct0 = new TDatime(2002,04,26,00,00,00);
+ TDatime ct = hitfile->GetCreationDate();
+
+ if(ct0->GetDate()>ct.GetDate()){
+ // For old files, must change SDD noise.
+ AliITSresponseSDD *resp1 = (AliITSresponseSDD*)ITS->DetType(1)->
+ GetResponseModel();
+ resp1 = new AliITSresponseSDD();
+ ITS->SetResponseModel(1,resp1);
+ cout << "Changed response class for SDD:" << endl;
+ resp1->Print();
+ } // end if
+
+ if(opt.Contains("Dubna")){
+ AliITSresponseSPDdubna *resp0 = new AliITSresponseSPDdubna();
+ if(ITS->DetType(0)->GetResponseModel() !=0){
+ delete ((AliITSresponse*)ITS->DetType(0)->GetResponseModel());
+ ITS->DetType(0)->ResponseModel(0);
+ } // end if
+ ITS->DetType(0)->ResponseModel(resp0);
+ AliITSsegmentationSPD *seg0 = (AliITSsegmentationSPD*)ITS->
+ DetType(0)->GetSegmentationModel();
+ AliITSsimulationSPDdubna *sim0 = new AliITSsimulationSPDdubna(seg0,
+ resp0);
+ if(ITS->DetType(0)->GetSimulationModel() !=0){
+ delete ((AliITSsimulation*)ITS->DetType(0)->GetSimulationModel());
+ ITS->DetType(0)->SimulationModel(0);
+ } // end if
+ ITS->DetType(0)->SimulationModel(sim0);
+ } // end if Dubna
+}
--- /dev/null
+// Neural performance evaluation
+//
+// before using this macro, you should have make run the macro
+// 'ITSstoreFindableTracks.C' in order to have a file named following
+// the rule to append '_fnd' to the name
+// (e.g. galice.root --> galice_fnd.root)
+// and have in it a tree with some tracks statistics useful for evaluation
+// of performance with the same event without loosing too much time...
+
+void AliITSNeuralCompleteEval
+(Int_t nsecs, const char *filename, Bool_t low = kFALSE,
+ Bool_t draw = kTRUE, const char *save = 0)
+{
+ Int_t N, M;
+ Float_t *edge = 0;
+ if (!low) {
+ N = 7;
+ edge = new Float_t[8];
+ edge[0] = 0.0;
+ edge[1] = 0.5;
+ edge[2] = 1.0;
+ edge[3] = 1.5;
+ edge[4] = 2.0;
+ edge[5] = 3.0;
+ edge[6] = 4.0;
+ edge[7] = 5.0;
+ }
+ else {
+ N = 5;
+ edge = new Float_t[6];
+ edge[0] = 0.0;
+ edge[1] = 0.2;
+ edge[2] = 0.4;
+ edge[3] = 0.6;
+ edge[4] = 0.8;
+ edge[5] = 1.0;
+ }
+
+ Float_t find[2] = {0.0, 0.0}, find1[2] = {0.0, 0.0};
+ Float_t good[2] = {0.0, 0.0}, good1[2] = {0.0, 0.0};
+ Float_t fake[2] = {0.0, 0.0}, fake1[2] = {0.0, 0.0};
+
+ // histos filled with neural results
+ TH1D *hgood[2], *hfake[2], *hfind[2];
+ TH1D *hgood1[2], *hfake1[2], *hfind1[2];
+ TH1D *hg[2], *hf[2];
+ if (draw) {
+ hgood[0] = new TH1D("hgood0", "Good found tracks", N, edge);
+ hfake[0] = new TH1D("hfake0", "Fake found tracks", N, edge);
+ hfind[0] = new TH1D("hfound0", "Findable tracks", N, edge);
+ hgood[1] = new TH1D("hgood1", "Good found tracks", N, edge);
+ hfake[1] = new TH1D("hfake1", "Fake found tracks", N, edge);
+ hfind[1] = new TH1D("hfound1", "Findable tracks", N, edge);
+
+ hgood[0]->Sumw2();
+ hfake[0]->Sumw2();
+ hfind[0]->Sumw2();
+ hgood[1]->Sumw2();
+ hfake[1]->Sumw2();
+ hfind[1]->Sumw2();
+
+ // histos for evaluating percentual efficiency
+ hg[0] = new TH1D("hg0", "Efficiency (%) for #geq 5 right pts.", N, edge);
+ hf[0] = new TH1D("hf0", "Fake probability (%) for #geq 5 right pts.", N, edge);
+ hg[1] = new TH1D("hg1", "Efficiency (%) for 6 right pts.", N, edge);
+ hf[1] = new TH1D("hf1", "Fake probability (%) for 6 right pts.", N, edge);
+ }
+
+ TFile *ffind;
+ TTree *tree;
+
+ if (low) {
+ ffind = new TFile(Form("%s.root", filename), "READ");
+ tree = (TTree*)ffind->Get("TreeF");
+ }
+ else {
+ ffind = new TFile(Form("%s_fnd.root", filename), "READ");
+ tree = (TTree*)ffind->Get("tree");
+ }
+
+ TFile *ftracks = new TFile(Form("%s_%d.root", filename, nsecs), "READ");
+
+
+ Double_t pt;
+ Int_t i, j, count, prim, *none = 0, div;
+ Int_t entries = tree->GetEntries(), label, min[] = {5,6};
+ tree->SetBranchAddress("pt", &pt);
+ tree->SetBranchAddress("count", &count);
+ tree->SetBranchAddress("prim", &prim);
+
+ AliITSneuralTrack *trk = 0;
+ div = low ? 50 : 500;
+ for(i = 1;;i++) {
+ trk = (AliITSneuralTrack*)ftracks->Get(Form("AliITSneuralTrack;%d", i));
+ if (i%div == 0) cout << "\rEvaluating found track " << i << flush;
+ if (!trk) break;
+ for (j = 0; j < 2; j++) {
+ label = trk->EvaluateTrack(0, min[j], none);
+ tree->GetEntry(abs(label));
+ if (count >= min[j]) {
+ if (label > 0) {
+ if (draw) hgood[j]->Fill(pt);
+ good[j]++;
+ if (pt >= 1.0) good1[j]++;
+ }
+ else {
+ if (draw) hfake[j]->Fill(pt);
+ fake[j]++;
+ if (pt >= 1.0) fake1[j]++;
+ }
+ }
+ }
+ }
+ cout << endl;
+
+ div = low ? 200 : 20000;
+
+ for (i = 0; i < entries; i++) {
+ if (i%div == 0) cout << "\rEvaluating findable track no. " << i << flush;
+ tree->GetEntry(i);
+ for (j = 0; j < 2; j++) {
+ if (count >= min[j]) {
+ find[j]++;
+ if (draw) hfind[j]->Fill(pt);
+ if (pt >= 1.0) find1[j]++;
+ }
+ }
+ }
+ cout << endl;
+ cout << hgood[0]->GetEntries() << " " << hgood[1]->GetEntries() << endl;
+ cout << hfake[0]->GetEntries() << " " << hfake[1]->GetEntries() << endl;
+ cout << hfind[0]->GetEntries() << " " << hfind[1]->GetEntries() << endl << endl;
+
+ if (draw) {
+ TCanvas *canv[2];
+ canv[0] = new TCanvas("c_0", "Tracking efficiency (soft)", 0, 0, 600, 500);
+ canv[1] = new TCanvas("c_1", "Tracking efficiency (hard)", 630, 0, 600, 500);
+
+ TLine *line1 = new TLine(1,100.0,edge[N],100.0); line1->SetLineStyle(4);
+ TLine *line2 = new TLine(1,90,edge[N],90); line2->SetLineStyle(4);
+
+ Bool_t good_drawn;
+ for (i = 0; i < 2; i++) {
+ canv[i]->cd();
+ good_drawn = kFALSE;
+ if (hgood[i]->GetEntries() > 0.0) {
+ good_drawn = kTRUE;
+ hg[i]->Divide(hgood[i], hfind[i], 100.0, 1.0);
+ hg[i]->SetMaximum(120);
+ hg[i]->SetMinimum(0);
+ hg[i]->SetMarkerStyle(21);
+ hg[i]->SetMarkerSize(1);
+ hg[i]->SetStats(kFALSE);
+ hg[i]->GetXaxis()->SetTitle("pt (GeV/c)");
+ hg[i]->Draw("PE1");
+ }
+ if (hfake[i]->GetEntries() > 0.0) {
+ hf[i]->Divide(hfake[i], hfind[i], 100.0, 1.0);
+ hf[i]->SetMaximum(120);
+ hf[i]->SetMinimum(0);
+ hf[i]->SetMarkerStyle(25);
+ hf[i]->SetMarkerSize(1);
+ hf[i]->SetStats(kFALSE);
+ if (good_drawn)
+ hf[i]->Draw("PE1SAME");
+ else
+ hf[i]->Draw("PE1");
+ }
+ line1->Draw("histosame");
+ line2->Draw("histosame");
+ canv[i]->Update();
+ }
+ canv[0]->SaveAs(Form("%s_soft.eps", filename));
+ canv[1]->SaveAs(Form("%s_hard.eps", filename));
+ cout << endl;
+ }
+
+ Float_t sgood[2] = {0.0, 0.0}, sgood1[2] = {0.0, 0.0};
+ Float_t sfake[2] = {0.0, 0.0}, sfake1[2] = {0.0, 0.0};
+ for (i = 0; i < 2; i++) {
+ sgood[i] = error(good[i], find[i]);
+ sgood1[i] = error(good1[i], find1[i]);
+ sfake[i] = error(fake[i], find[i]);
+ sfake1[i] = error(fake1[i], find1[i]);
+
+ good[i] = good[i] * 100.0 / find[i];
+ fake[i] = fake[i] * 100.0 / find[i];
+ good1[i] = good1[i] * 100.0 / find1[i];
+ fake1[i] = fake1[i] * 100.0 / find1[i];
+ }
+
+ if (save) {
+ fstream data(save, ios::app);
+ data.setf(ios::fixed);
+ data.precision(1);
+ data << good1[0] << " " << fake1[0] << " ";
+ data << good1[1] << " " << fake1[1] << endl;
+ data.close();
+ }
+ else {
+ cout.setf(ios::fixed);
+ cout.precision(1);
+ cout << "*****************************************" << endl;
+ cout << "* Tracks with at least 5 correct points *" << endl;
+ cout << "*****************************************" << endl;
+ cout << "(all particles)" << endl;
+ cout << "Efficiency: " << good[0] << " +/- " << sgood[0] << "%" << endl;
+ cout << "Fake prob.: " << fake[0] << " +/- " << sfake[0] << "%" << endl;
+ if (!low) {
+ cout << "(pt >= 1 GeV/c)" << endl;
+ cout << "Efficiency: " << good1[0] << " +/- " << sgood1[0] << "%" << endl;
+ cout << "Fake prob.: " << fake1[0] << " +/- " << sfake1[0] << "%" << endl;
+ }
+ cout << endl;
+ cout << "************************************" << endl;
+ cout << "* Tracks with all 6 correct points *" << endl;
+ cout << "************************************" << endl;
+ cout << "(all particles)" << endl;
+ cout << "Efficiency: " << good[1] << " +/- " << sgood[1] << "%" << endl;
+ cout << "Fake prob.: " << fake[1] << " +/- " << sfake[1] << "%" << endl;
+ if (!low) {
+ cout << "(pt >= 1 GeV/c)" << endl;
+ cout << "Efficiency: " << good1[1] << " +/- " << sgood1[1] << "%" << endl;
+ cout << "Fake prob.: " << fake1[1] << " +/- " << sfake1[1] << "%" << endl;
+ }
+ cout << endl;
+ }
+}
+
+Double_t error(Double_t x, Double_t y) {
+ Double_t radq = x + x * x / y;
+ return 100.0 * sqrt(radq) / y;
+}
--- /dev/null
+// ARGUMENTS:
+// 1. number of azymuthal sectors (it's better not to go under 8 or over 40)
+// 2. the ROOT file to read (WITHOUT exstension)
+// 3. event number
+// 4. if specified a string, a fstream named like the argument is opened and
+// the elapsed CPU time is stored (not useful)
+
+// the macro will save a file named, for example "galice_<nsecs>.root"
+// containing may AliITSneuralTrack objects
+
+void AliITSNeuralTracking
+(Int_t nsecs = 12, const char* rfile = "galice", Int_t event = 0, const char* save = 0)
+{
+ TStopwatch timer;
+ Double_t CONVERT = TMath::Pi() / 180.0;
+ const char* wfile = Form("%s_%d.root", rfile, nsecs);
+ cout << "Reading file " << rfile << ".root and saving in " << wfile << endl;
+
+// ==================================
+// ==== CURVATURE CUT DEFINITION ====
+// ==================================
+
+ // These values define the curvature cuts for all steps
+ // within a sector.
+ // For a greater clarity, the cuts are indicated in units
+ // of transverse momentum (GeV/c) but these value have no
+ // exact physical meaning, but are useful to understand
+ // well what means a choice in the value of a certain
+ // curvature constraint
+ // NOTE: becareful to make sure that the 'ncuts' variable
+ // have the same value of the dimension of the allocated arrays
+
+ Int_t ncuts;
+ Double_t *p, *cut;
+
+ ncuts = 5;
+ p = new Double_t[5];
+ cut = new Double_t[5];
+ p[0] = 2.0;
+ p[1] = 1.0;
+ p[2] = 0.7;
+ p[3] = 0.5;
+ p[4] = 0.3;
+
+ for (Int_t i = 0; i < ncuts; i++) cut[i] = 0.003 * 0.2 / p[i];
+
+
+// ==========================
+// ==== OTHER PARAMETERS ====
+// ==========================
+
+ Bool_t flag = kFALSE; // for now, don't change this line, please...
+
+ Double_t diff = 0.02; // helicoidal cut
+ Double_t dtheta = 1.0; // delta-theta cut
+ Double_t temp = 1.0; // temperature parameter
+ Double_t var = 0.0001; // stabilization threshold
+
+ Double_t exp = 7.0; // straight-line excitator
+ Double_t gtoc = 3.0; // gain/cost contribution ratio
+
+ Double_t min = 0.4; // minimum in random activation initialization
+ Double_t max = 0.6; // maximum in random activation initialization
+
+
+// =========================
+// ==== NEURAL TRACKING ====
+// =========================
+
+ AliITSneuralTracker *ANN = new AliITSneuralTracker(nsecs, ncuts, cut, CONVERT*dtheta);
+
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(Form("%s.root", rfile));
+ if (!file) file = new TFile(Form("%s.root", rfile),"UPDATE");
+
+ //Double_t Xv = -0.001097;
+ //Double_t Yv = -0.00347647;
+ //Double_t Zv = 0.000631345;
+ //ANN->SetVertex(Xv, Yv, Zv);
+ // You should find the vertex with VertexMacro.C
+ // and then put by hand the found values with
+ // the above method.
+
+ Int_t points = ANN->ReadFile(Form("%s.root", rfile), event);
+
+ ANN->SetTemperature(temp);
+ ANN->SetVariationLimit(var);
+ ANN->SetGainToCostRatio(gtoc);
+ ANN->SetExponent(exp);
+ ANN->SetInitLimits(min, max);
+ ANN->SetDiff(diff);
+
+ cout << points << " points found " << endl << endl;
+
+ TStopwatch timer;
+ timer.Start();
+
+ ANN->Go(wfile, flag);
+
+ timer.Stop();
+ cout << endl;
+ timer.Print();
+
+ if (save) {
+ fstream ksave(save, ios::app);
+ ksave << nsecs << " " << timer->CpuTime() << endl;
+ ksave.close();
+ }
+
+// delete gAlice;
+// gAlice = 0;
+}
--- /dev/null
+/******************************************************************************
+
+ "AliITSOccupancy.C"
+
+ this macro calculates the mean occupancy of each ITS layer, counting the
+ "fired" digits of each module, and making two overall calculations:
+ 1) the calculus of the mean overall layer occupancy (as the ratio
+ between the total number of active digits and the total number of
+ channels in the module;
+ 2) a distribution of the occupancies as a funcion of z, to obtain some
+ kind of most significand data for this value, along the z axis
+
+ The macro creates also the GIF and the EPS of the TCanvas where the
+ histograms are plotted. Muliple input files are supported (see the
+ input argument list below).
+
+
+ It is also possible (and recommended) to compile this macro,
+ to execute it faster. To do this, it is necessary to:
+
+ 1) type, at the root prompt, the instruction
+ gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g")
+ to adjoin the ALICE include libraries to the main include directory
+ 2) load the function via the ".L AliITSOccupancy.C++" statement
+ (only the first time, because the first time the macro must be compiled).
+ From the second time you use the macro, you must only load it via the
+ ".L AliITSOccupancy.C+" instruction (note the number of '+' signs in
+ each case
+ 3) call the function with only its name, e.g.: AliITSOccupancy()
+
+ Author: Alberto Pulvirenti
+******************************************************************************/
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+ #include <iostream.h>
+ #include <TROOT.h>
+ #include <TArrayI.h>
+ #include <TCanvas.h>
+ #include <TClassTable.h>
+ #include <TClonesArray.h>
+ #include <TFile.h>
+ #include <TObject.h>
+ #include <TObjArray.h>
+ #include <TTree.h>
+ #include <TMath.h>
+ #include <TString.h>
+ #include <TH1.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>
+#endif
+
+TFile* AccessFile(TString inFile="galice.root", TString acctype="R");
+
+Int_t AliITSOccupancy(TString FileHits="galice.root", TString FileDigits="galice.root", Int_t evNum = 0) {
+
+ // Open the Root input file containing the AliRun data and
+ // the Root output data file
+ TFile *froot = AccessFile(FileHits);
+ froot->ls();
+
+ if(!(FileDigits.Data() == FileHits.Data())){
+ gAlice->SetTreeDFileName(FileDigits);
+ }
+ gAlice->GetEvent(evNum);
+ // Initialize the ITS and its geometry
+ AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
+ AliITSgeom *gm = ITS->GetITSgeom();
+
+ // Variables and objects definition and setting
+ Float_t zmax[] = {20.0,20.0,25.0,35.0,49.5,54.0}; // z edges for TH1F (cm)
+ Int_t nbins[] = {40,40,50,70,99,108}; // bins number for TH1Fs
+
+ // Histos for plotting occupancy distributions
+ TH1F *above[6], *below[6];
+
+ for (Int_t lay = 0; lay < 6; lay++) {
+ Int_t nlads = gm->GetNladders(lay+1);
+ Int_t ndets = gm->GetNdetectors(lay+1);
+ Int_t dtype = lay / 2;
+ Int_t minID = gm->GetModuleIndex(lay+1, 1, 1);
+ Int_t maxID = gm->GetModuleIndex(lay+1, nlads, ndets);
+ Text_t ffname[20];
+ sprintf(ffname, "h_%d", lay+1);
+ below[lay] = new TH1F(ffname, "Total z distribution of digits", nbins[lay], -zmax[lay], zmax[lay]);
+ cout << "Evaluating total channels number of layer " << lay+1 << endl;
+ for (Int_t I = minID; I <= maxID; I++) {
+ AliITSDetType *det = ITS->DetType(dtype);
+ AliITSsegmentation *seg = det->GetSegmentationModel();
+ Int_t NX = seg->Npx();
+ if(lay == 2 || lay == 3) NX = 192;
+ Int_t NZ = seg->Npz();
+ // cout << "ID: " << I << ", Layer: " << lay+1 << ", NX = " << NX << ", NZ = " << NZ << endl;
+ for (Int_t ix = 0; ix <= NX; ix++) {
+ for (Int_t iz = 0; iz <= NZ; iz++) {
+ Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
+ seg->DetToLocal(ix, iz, lx[0], lx[2]);
+ gm->LtoG(I, lx, gx);
+ below[lay]->Fill(gx[2]);
+ }
+ }
+ }
+ sprintf(ffname, "H_%d", lay+1);
+ above[lay] = new TH1F(ffname, "histo", nbins[lay], -zmax[lay], zmax[lay]);
+ }
+
+ // Counting the hits, digits and recpoints contained in the ITS
+ TTree *TD = gAlice->TreeD();
+ ITS->ResetDigits();
+ Float_t mean[6];
+ Float_t average[6];
+
+ for (Int_t L = 0; L < 6; L++) {
+
+ cout << "Layer " << L + 1 << ": " << flush;
+
+ // To avoid two nested loops, are calculated
+ // the ID of the first and last module of the L
+ // (remember the L goest from 0 to 5 (not from 1 to 6)
+ Int_t first, last;
+ first = gm->GetModuleIndex(L + 1, 1, 1);
+ last = gm->GetModuleIndex(L + 1, gm->GetNladders(L + 1), gm->GetNdetectors(L + 1));
+
+ // Start loop on modules
+ for (Int_t ID = first; ID <= last; ID++) {
+ Int_t dtype = L / 2;
+ AliITSDetType *det = ITS->DetType(dtype);
+ AliITSsegmentation *seg = det->GetSegmentationModel();
+ if (dtype == 2) seg->SetLayer(L+1);
+
+ TD->GetEvent(ID);
+ TClonesArray *digits_array = ITS->DigitsAddress(dtype);
+ Int_t digits_num = digits_array->GetEntries();
+ // Get the coordinates of the module
+ for (Int_t j = 0; j < digits_num; j++) {
+ Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
+ AliITSdigit *digit = (AliITSdigit*)digits_array->UncheckedAt(j);
+ Int_t iz=digit->fCoord1; // cell number z
+ Int_t ix=digit->fCoord2; // cell number x
+ // Get local coordinates of the element (microns)
+ seg->DetToLocal(ix, iz, lx[0], lx[2]);
+ gm->LtoG(ID, lx, gx);
+ above[L]->Fill(gx[2]);
+ }
+ }
+
+ Float_t occupied = above[L]->GetEntries();
+ Float_t total = below[L]->GetEntries();
+ cout << "Entries: " << occupied << ", " << total << endl;
+ average[L] = 100.*occupied/total;
+ above[L]->Divide(above[L], below[L], 100.0, 1.0);
+ mean[L] = above[L]->GetSumOfWeights() / above[L]->GetNbinsX();
+ cout.setf(ios::fixed);
+ cout.precision(2);
+ cout << " digits occupancy = " << mean[L] << "%" << endl;
+ cout << " average digits occupancy = " << average[L] << "%" << endl;
+ }
+
+ TCanvas *view = new TCanvas("view", "Digits occupancy distributions", 600, 900);
+ view->Divide(2, 3);
+
+ for (Int_t I = 1; I < 7; I++) {
+ view->cd(I);
+ Text_t title[50];
+ sprintf(title, "Layer %d: %4.2f%c", I, mean[I-1], '%');
+ title[6] = (Text_t)I + '0';
+ above[I-1]->SetTitle(title);
+ above[I-1]->SetStats(kFALSE);
+ above[I-1]->SetXTitle("z (cm)");
+ above[I-1]->SetYTitle("%");
+ above[I-1]->Draw();
+ view->Update();
+ }
+
+ view->SaveAs("AliITSOccupancy_digit.gif");
+ view->SaveAs("AliITSOccupancy_digit.eps");
+ TFile *fOc = new TFile("AliITSOccupancy.root","recreate");
+ fOc->cd();
+ for (Int_t I = 0; I < 6; I++) above[I]->Write();
+ fOc->Close();
+
+ return 1;
+}
+
+//-------------------------------------------------------------------
+TFile * AccessFile(TString FileName, TString acctype){
+
+ // Function used to open the input file and fetch the AliRun object
+
+ TFile *retfil = 0;
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(FileName);
+ if (file) {file->Close(); delete file; file = 0;}
+ if(acctype.Contains("U")){
+ file = new TFile(FileName,"update");
+ }
+ if(acctype.Contains("N") && !file){
+ file = new TFile(FileName,"recreate");
+ }
+ if(!file) file = new TFile(FileName); // default readonly
+ if (!file->IsOpen()) {
+ cerr<<"Can't open "<<FileName<<" !" << endl;
+ return retfil;
+ }
+
+ // Get AliRun object from file or return if not on file
+ if (gAlice) {delete gAlice; gAlice = 0;}
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (!gAlice) {
+ cerr << "AliRun object not found on file"<< endl;
+ return retfil;
+ }
+ return file;
+}
--- /dev/null
+void AliITSPlotGeom(const char *filename="galice.root"){
+/////////////////////////////////////////////////////////////////////////
+// This macro displays part of the ITS sensitive volume as defined
+// in the root file.
+//
+// Root > .L AliITSPlotGeom.C //this loads the macro in memory
+// Root > AliITSPlotGeom(); //by default process first event
+// Root > AliITSPlotGeom("galice2.root"); //process third event from
+// galice2.root file.
+//Begin_Html
+/*
+<img src="picts/AliITSplotgeom.gif">
+*/
+//End_Html
+/////////////////////////////////////////////////////////////////////////
+ if(gAlice){
+ delete gAlice;
+ gAlice=0;
+ }else{
+ // Dynamically link some shared libs
+ if(gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ } // end if
+ } // end if gAlice
+ // Connect the Root Galice file containing Geometry, Kine and Hits
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
+ if(!file) file = new TFile(filename);
+
+ // 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");
+ } // end if !gAlice
+
+
+// Get pointers to the ITS Alice detectors and Hits containers
+ AliITS *ITS = (AliITS*) gAlice->GetDetector("ITS");
+ if(ITS==0){
+ cout << "ITS not found. Exiting." << endl;
+ return;
+ } // end if ITS==0
+ AliITSgeom *itsgeom = ITS->GetITSgeom();
+ if(itsgeom==0){
+ cout << "ITSgeom not defined. Exiting." << endl;
+ return;
+ } // end if itsgeom==0
+ Int_t is,ie,in;
+ is = itsgeom->GetStartSPD();
+ ie = itsgeom->GetLastSPD()+1;
+ in = ie-is;
+
+ // Transformations.
+// Float_t *x = new Float_t[in];
+// Float_t *y = new Float_t[in];
+// Float_t *z = new Float_t[in];
+ Float_t x,y,z;
+ TRotMatrix **r = new TRotMatrix*[in];
+ Int_t i,j;
+ Double_t m[9];
+ char name[10],title[20],name2[10],title2[20];
+// cout << in << endl;
+ for(i=0;i<in;i++){
+ itsgeom->GetRotMatrix(i+is,m);
+ sprintf(name,"ROT%d",i+is);
+ sprintf(title,"ROT%d",i+is);
+ r[i] = new TRotMatrix(name,title,m);
+// cout << name << title << endl;
+// itsgeom->GetTrans(i+is,x[i],y[i],z[i]);
+// cout << i << " " << x[i] << " " << y[i] << " " << z[i] <<endl;
+ } // end for i
+
+ // Sensitive volumes.
+ if(itsgeom->IsShapeDefined(0)){
+ TBRIK *spds = (TShape*) (itsgeom->GetShape(0));
+ }else{
+ TBRIK *spds = new TBRIK("SPD","SPD","void",0.64,0.015,3.48);
+ } // end if
+// cout << spds << endl;
+ if(itsgeom->IsShapeDefined(1)){
+ TBRIK *sdds = (TShape*) (itsgeom->GetShape(240));
+ }else{
+ TBRIK *sdds = new TBRIK("SDD","SDD","void",3.50850,0.01499,3.76320);
+ } // end if
+// cout << sdds << endl;
+ if(itsgeom->IsShapeDefined(2)){
+ TBRIK *ssds = (TShape*) (itsgeom->GetShape(1000));
+ }else{
+ TBRIK *ssds = new TBRIK("SSD","SSD","void",3.65,0.015,2.00);
+ } // end if
+// cout << ssds << endl;
+
+ // Set up display.
+
+ TCanvas *c1 = new TCanvas("c1","ITS geometry",10,10,700,700);
+// TPad *p1 = new TPad("p1","p1",0.01,0.01,0.99,0.99,46,3,1);
+// p1->Draw();
+// p1->cd();
+// TView *view = new TView(1);
+// view->SetRange(-5,-5,-5,5,5,5);
+ TShape *mother = new TBRIK("Mother","Mother Volume","void",10,10,10);
+// TShape *mother = new TBRIK("Mother","Mother Volume","void",30,30,30);
+// TShape *mother = new TBRIK("Mother","Mother Volume","void",50,50,50);
+ TNode *node = new TNode("node","Mother Node",mother);
+ node->cd();
+
+ // Set up nodes
+ TNode **itsn = new TNode*[in];
+ TPolyLine3D **pl = new TPolyLine3D*[in];
+/* Double_t p[5*3],axis[5*3]={1.,0.,0.,
+ 0.,0.,0.,
+ 0.,1.,0.,
+ 0.,0.,0.,
+ 0.,0.,1.}; */
+ Double_t p[19*3],axis[19*3]={-0.25,0.,1.25, // 1
+ +0.00,0.,1.25, // 2
+ -0.25,0.,1.00, // 3
+ +0.00,0.,1.00, // 4
+ +0.00,0.,0.00, // 5
+ +1.00,0.,0.00, // 6
+ +1.25,0.,0.25, // 7
+ +1.125,0.,0.125, // 8
+ +1.00,0.,0.25, // 9
+ +1.125,0.,0.125, // 10
+ +1.25,0.,0.00, // 11
+ +1.125,0.,0.125, // 12
+ +1.00,0.,0.00, // 13
+ +0.00,0.,0.00, // 14
+ +0.00,1.,0.00, // 15
+ +0.00,1.,0.125, // 16
+ +0.25,1.,0.25, // 17
+ +0.00,1.,0.125, // 18
+ -0.25,1.,0.25, // 19
+ };
+ for(i=0;i<in;i++){
+ itsgeom->GetTrans(i+is,x,y,z);
+/* switch (itsgeom->GetGeomMatrix(i+is)->GetDetectorIndex())
+ case 0:
+ sprintf(name,"SPD%d",i+is);
+ sprintf(title,"SPD%d",i+is);
+ sprintf(name2,"BSPD%d",i+is);
+ sprintf(title2,"BSPD%d",i+is);
+ itsn[i] = new TNode(name,title,new TBRIK(name2,title2,"void",
+ spds->GetDx(),spds->GetDy(),spds->GetDz()),x,y,z,r[i]);
+ break;
+ case 1:
+ sprintf(name,"SDD%d",i+is);
+ sprintf(title,"SDD%d",i+is);
+ sprintf(name2,"BSDD%d",i+is);
+ sprintf(title2,"BSDD%d",i+is);
+ itsn[i] = new TNode(name,title,new TBRIK(name2,title2,"void",
+ sdds->GetDx(),sdds->GetDy(),sdds->GetDz()),x,y,z,r[i]);
+ break;
+ case 2:
+ sprintf(name,"SSD%d",i+is);
+ sprintf(title,"SSD%d",i+is);
+ sprintf(name2,"BSSD%d",i+is);
+ sprintf(title2,"BSSD%d",i+is);
+ itsn[i] = new TNode(name,title,new TBRIK(name2,title2,"void",
+ ssds->GetDx(),ssds->GetDy(),ssds->GetDz()),x,y,z,r[i]);
+ break; */
+ for(j=0;j<19;j++) itsgeom->LtoG(i+is,(Double_t*)&axis[3*j],
+ (Double_t*)&p[3*j]);
+ pl[i] = new TPolyLine3D(19,p);
+ } // end for i
+
+ // display it
+ node->cd();
+ node->Draw();
+// for(i=0;i<in;i++) itsn[i]->Draw();
+// node->Draw();
+ for(i=0;i<in;i++) pl[i]->Draw();
+ c1->Update();
+
+ // clean up
+// delete[] x;
+// delete[] y;
+// delete[] z;
+// for(i=0;i<itsgeom->GetIndexMax();i++) delete[] r[i];
+// delete[] r;
+}
+
--- /dev/null
+#include <iostream.h>
+#include <fstream.h>
+
+void AliITSPlotTracksV1()
+{
+ TFile *file = new TFile("AliITSComparisonV1.root");
+ TTree *tree = (TTree*)file->Get("Eval");
+
+ if (!tree) {
+ cerr << "Evaluation tree not found!" << endl;
+ return;
+ }
+
+ // histogram definition - I (efficiency)
+ TH1D *hFindables = new TH1D("hFindables", "Findable tracks", 10, 0, 2);
+ TH1D *hGood = new TH1D("hGood", "Good found tracks", 10, 0, 2);
+ TH1D *hFake = new TH1D("hFake", "Fake found tracks", 10, 0, 2);
+ TH1D *hRatioG = new TH1D("hRatioG", "", 10, 0, 2);
+ TH1D *hRatioF = new TH1D("hRatioF", "", 10, 0, 2);
+ hGood->Sumw2();
+ hFake->Sumw2();
+ hFindables->Sumw2();
+ hRatioG->SetLineColor(kBlue);
+ hRatioG->SetLineWidth(2);
+ hRatioF->SetLineColor(kRed);
+ hRatioF->SetLineWidth(2);
+
+ // histograms definition - II (resolution)
+ TH1D *hPhi = new TH1D("hPhi", "#Phi resolution", 50, -15., 15.);
+ TH1D *hLambda = new TH1D("hLambda", "#lambda resolution", 50, -15., 15.);
+ TH1D *hPt = new TH1D("hPt", "Relative Pt resolution", 40, -10., 10.);
+ TH1D *hDtot = new TH1D("hDtot", "Total impact parameter distribution", 100, 0., 2000.);
+ TH1D *hDr = new TH1D("hDr", "Transv. impact parameter distribution", 50, -1000., 1000.);
+ TH1D *hDz = new TH1D("hDz", "Long. impact parameter distribution", 50, -1000., 1000.);
+ hPhi->SetFillColor(4);
+ hLambda->SetFillColor(4);
+ hPt->SetFillColor(2);
+ hDtot->SetFillColor(6);
+ hDr->SetFillColor(kGreen);
+ hDz->SetFillColor(kBlue);
+
+ // Evaluation tree settings
+ Int_t labITS, labTPC, signC;
+ Int_t i, j, tot = (Int_t)tree->GetEntries();
+ Double_t difpt, diflambda, difphi, Dz, Dr, Dtot, ptg;
+ tree->SetBranchAddress("labITS" , &labITS );
+ tree->SetBranchAddress("labTPC" , &labTPC );
+ tree->SetBranchAddress("difpt" , &difpt );
+ tree->SetBranchAddress("diflambda", &diflambda);
+ tree->SetBranchAddress("difphi" , &difphi );
+ tree->SetBranchAddress("Dz" , &Dz );
+ tree->SetBranchAddress("Dr" , &Dr );
+ tree->SetBranchAddress("Dtot" , &Dtot );
+ tree->SetBranchAddress("ptg" , &ptg );
+ tree->SetBranchAddress("signC" , &signC );
+
+ // Filling the histogram of findable tracks (w.r. to momentum)
+ for(i = 0; i < tot; i++) {
+ tree->GetEntry(i);
+ hFindables->Fill(ptg);
+ }
+
+ // Filling the evaluation histograms
+ Int_t neglabs = 0;
+ for(i = 0; i < tot; i++) {
+ tree->GetEntry(i);
+// if(signC<0) continue;
+ if (labITS < 0) neglabs++;
+ if (labITS >= 0) {
+ hGood->Fill(ptg);
+ hPt->Fill(difpt);
+ hLambda->Fill(diflambda);
+ hPhi->Fill(difphi);
+ hDtot->Fill(Dtot);
+ hDr->Fill(Dr);
+ hDz->Fill(Dz);
+ }
+ else {
+ hFake->Fill(ptg);
+ neglabs++;
+ }
+ }
+
+ // Drawing
+ cerr << "Findable tracks : " << hFindables->GetEntries() << endl;
+ cerr << "Good found tracks : " << hGood->GetEntries() << endl;
+ cerr << "Fake found tracks : " << hFake->GetEntries() << endl;
+
+ gStyle->SetOptStat(111110);
+ gStyle->SetOptFit(1);
+
+ TCanvas *c1 = new TCanvas("c1","Parameter resolutions",0,0,700,700);
+ c1->Divide(2, 2, 0.001, 0.001);
+ c1->cd(1); hPhi->SetXTitle("(mrad)"); hPhi->Draw(); hPhi->Fit("gaus", "Q");
+ c1->cd(2); hLambda->SetXTitle("(mrad)"); hLambda->Draw(); hLambda->Fit("gaus", "Q");
+ c1->cd(3); hPt->SetXTitle("(%)"); hPt->Draw(); hPt->Fit("gaus", "Q");
+ c1->cd(4); hDtot->SetXTitle("(micron)"); hDtot->Draw();
+ c1->Update();
+
+ TCanvas *c2 = new TCanvas("c2","Impact parameters resolution",100,100,700,400);
+ c2->Divide(2, 1, 0.001, 0.001);
+ c2->cd(1); hDr->SetXTitle("(micron)"); hDr->Draw(); hDr->Fit("gaus", "Q");
+ c2->cd(2); hDz->SetXTitle("(micron)"); hDz->Draw(); hDz->Fit("gaus", "Q");
+ c2->Update();
+
+ TCanvas *c3 = new TCanvas("c3","Momentum distributions",200,200,800,500);
+ c3->Divide(3, 1, 0.001, 0.001);
+ c3->cd(1); hGood->Draw();
+ c3->cd(2); hFake->Draw();
+ c3->cd(3); hFindables->Draw();
+ c3->Update();
+
+ TCanvas *c4 = new TCanvas("c4","Tracking efficiency",300,300,800,500);
+ hRatioG->Divide(hGood, hFindables, 1., 1.);
+ hRatioF->Divide(hFake, hFindables, 1., 1.);
+ hRatioG->SetMaximum(1.4);
+ hRatioG->SetYTitle("Tracking efficiency");
+ hRatioG->SetXTitle("Pt (GeV/c)");
+ hRatioG->Draw(); // to not plot the erro bars hg->Draw("histo");
+ hRatioF->SetFillColor(1);
+ hRatioF->SetFillStyle(3013);
+ hRatioF->SetLineColor(2);
+ hRatioF->SetLineWidth(2);
+ hRatioF->Draw("same"); // to not plot the error bars hRatioF->Draw("histosame");
+ // a line to mark the best efficiency
+ TLine *line1 = new TLine(0,1.0,2,1.0); line1->SetLineStyle(4);
+ line1->Draw("same");
+ TLine *line2 = new TLine(0,0.9,2,0.9); line2->SetLineStyle(4);
+ line2->Draw("histosame");
+ // a text explanation
+ TText *text = new TText(0.461176,0.248448,"Fake tracks");
+ text->SetTextSize(0.05);
+ text->Draw();
+ text = new TText(0.453919,1.11408,"Good tracks");
+ text->SetTextSize(0.05);
+ text->Draw();
+ c4->Update();
+
+ cout << "neglabs = " << neglabs << endl; // provvisoria
+}
--- /dev/null
+#include "iostream.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TClonesArray.h"
+/*
+#include "$(ALICE_ROOT)/STEER/AliRun.h"
+#include "$(ALICE_ROOT)/ITS/AliITS.h"
+#include "$(ALICE_ROOT)/ITS/AliITSgeom.h"
+#include "$(ALICE_ROOT)/ITS/AliITSdigit.h"
+*/
+void AliITSPrintDigits(TString rfn="galice.root",Int_t mod=-1,
+ Int_t evnt=-1){
+ // Macro to print out the recpoints for all or a specific module
+
+ // Dynamically link some shared libs
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ } // end if
+ gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSstandard.C");
+
+ TFile *rf=0;
+ rf = AccessFile(rfn,"R"); // Set up to read in Data
+ AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
+ if(!ITS){
+ cout << "Error: no ITS found. Aborting"<<endl;
+ return;
+ } // end if !ITS
+
+ Int_t evNumber1 = 0;
+ Int_t evNumber2 = gAlice->GetEventsPerRun();
+ if(evnt>=0){
+ evNumber1 = evnt;
+ evNumber2 = evnt+1;
+ } // end if evnt>=0
+ Int_t mod1 = 0;
+ Int_t mod2 = ITS->GetITSgeom()->GetIndexMax();
+ if(mod>=0){
+ mod1 = mod;
+ mod2 = mode+1;
+ } // end if mod>=0
+ TClonesArray *pda = 0;
+ AliITSdigitSPD *dp0 = 0;
+ AliITSdigitSDD *dp1 = 0;
+ AliITSdigitSSD *dp2 = 0;
+
+ Int_t event,m,i,i2,id;
+ for(event = evNumber1; event < evNumber2; event++){
+ gAlice->GetEvent(event);
+ for(m=mod1;m<mod2;m++){
+ id = ITS->GetITSgeom()->GetModuleType(m);
+ dpa = ITS->DigitsAddress(id);
+ ITS->ResetDigits();
+ gAlice->TreeD()->GetEvent(m);
+ i2 = dpa->GetEntriesFast();
+ switch (id) {
+ case 0:
+ cout << "Event=" << event << " module=" << m <<
+ " Number of SPD Digits=" << i2 <<endl;
+ for(i=0;i<i2;i++){
+ dp0 = (AliITSdigitSPD*)(dpa->At(i));
+ cout << i << " ";
+ dp0->Print((ostream*)cout);
+ cout << endl;
+ } // end for i
+ break;
+ case 1:
+ cout << "Event=" << event << " module=" << m <<
+ " Number of SDD Digits=" << i2 <<endl;
+ for(i=0;i<i2;i++){
+ dp1 = (AliITSdigitSDD*)(dpa->At(i));
+ cout << i << " ";
+ dp1->Print((ostream*)cout);
+ cout << endl;
+ } // end for i
+ break;
+ case 2:
+ cout << "Event=" << event << " module=" << m <<
+ " Number of SSD Digits=" << i2 <<endl;
+ for(i=0;i<i2;i++){
+ dp2 = (AliITSdigitSSD*)(dpa->At(i));
+ cout << i << " ";
+ dp2->Print((ostream*)cout);
+ cout << endl;
+ } // end for i
+ break;
+ default:
+ break;
+ } // end switch
+ } // end for m
+ } // end for event
+}
--- /dev/null
+#include "iostream.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TClonesArray.h"
+/*
+#include "$(ALICE_ROOT)/STEER/AliRun.h"
+#include "$(ALICE_ROOT)/ITS/AliITS.h"
+#include "$(ALICE_ROOT)/ITS/AliITSgeom.h"
+#include "$(ALICE_ROOT)/ITS/AliITSpList.h"
+*/
+void AliITSPrintSDigits(TString rfn="galice.root",Int_t mod=-1,
+ Int_t evnt=-1){
+ // Macro to print out the recpoints for all or a specific module
+
+ // Dynamically link some shared libs
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ } // end if
+ gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSstandard.C");
+
+ TFile *rf=0;
+ rf = AccessFile(rfn,"R"); // Set up to read in Data
+ AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
+ if(!ITS){
+ cout << "Error: no ITS found. Aborting"<<endl;
+ return;
+ } // end if !ITS
+
+ Int_t evNumber1 = 0;
+ Int_t evNumber2 = gAlice->GetEventsPerRun();
+ if(evnt>=0){
+ evNumber1 = evnt;
+ evNumber2 = evnt+1;
+ } // end if evnt>=0
+ Int_t mod1 = 0;
+ Int_t mod2 = ITS->GetITSgeom()->GetIndexMax();
+ if(mod>=0){
+ mod1 = mod;
+ mod2 = mode+1;
+ } // end if mod>=0
+ TClonesArray *sdpa = ITS->GetSDigits();
+ AliITSpListItem *sdp = 0;
+
+ Int_t event,m,i,i2;
+ for(event = evNumber1; event < evNumber2; event++){
+ gAlice->GetEvent(event);
+ for(m=mod1;m<mod2;m++){
+ ITS->ResetSDigits();
+ gAlice->TreeS()->GetEvent(m);
+ i2 = sdpa->GetEntriesFast();
+ cout << "Event=" << event << " module=" << m <<
+ " Number of SDigits=" << i2 <<endl;
+ for(i=0;i<i2;i++){
+ sdp = (AliITSpListItem*)(sdpa->At(i));
+ cout << i << " ";
+ sdp->Print((ostream*)cout);
+ cout << endl;
+ } // end for i
+ } // end for m
+ } // end for event
+
+}
--- /dev/null
+#ifndef __CINT__
+ #include <Riostream.h>
+ #include "AliITSgeom.h"
+ #include "AliITStrackerV2.h"
+
+ #include "TFile.h"
+ #include "TStopwatch.h"
+#endif
+
+Int_t AliITSPropagateBackV2(Int_t nev=1) {
+ cerr<<"Propagating tracks back through the ITS...\n";
+
+ TFile *in=TFile::Open("AliITStracksV2.root");
+ if (!in->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 1;}
+
+ TFile *out=TFile::Open("AliTPCtracks.root","update");
+ if (!out->IsOpen()) {
+ cerr<<"Can't open AliTPCtracks.root !\n"; return 2;
+ }
+ TFile *file=TFile::Open("AliITSclustersV2.root");
+ if (!file->IsOpen()) {
+ cerr<<"Can't open AliITSclustersV2.root !\n";return 3;
+ }
+ AliITSgeom *geom=(AliITSgeom*)file->Get("AliITSgeom");
+
+ Int_t rc=0;
+ TStopwatch timer;
+ AliITStrackerV2 tracker(geom);
+ for (Int_t i=0; i<nev; i++) {
+ cerr<<"Processing event number : "<<i<<endl;
+ tracker.SetEventNumber(i);
+ rc=tracker.PropagateBack(in,out);
+ }
+ timer.Stop(); timer.Print();
+
+ file->Close();
+ in->Close();
+ out->Close();
+
+ return rc;
+}
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include<iostream.h>
+#include<TROOT.h>
+#include<TArrayI.h>
+#include<TBranch.h>
+#include<TCanvas.h>
+#include<TClassTable.h>
+#include<TClonesArray.h>
+#include<TFile.h>
+#include<TH1.h>
+#include<TH2.h>
+#include<TLegend.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>
+#endif
+
+Int_t GetModuleHits (TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St);
+Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z, Int_t isfastpoints);
+Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z);
+
+Int_t AliITSReadPlotData( Int_t evNum = 0, char *filename = "galice.root", TString FileDigits="galice.root", TString FileRec="galice.root", Int_t isfastpoints = 0) {
+
+ /*********************************************************************
+ * *
+ * Macro used to read hits, digits and recpoints for a module *
+ * It draws 3 TH2Fs where stores the 2-D coordinates of these *
+ * data for a specified module (for which the macro asks when *
+ * it starts, with some checks to avoid wrong detector coords. *
+ * *
+ * Only a little 'experimental' warning: *
+ * with all the tests I have made, I wasn't able to plot the *
+ * digits fo the 5th layer... *
+ * I skipped this problem with an alert to beware th user, while *
+ * in this situation the macro plots only recpoints and hits *
+ * *
+ * Author: Alberto Pulvirenti *
+ * *
+ *********************************************************************/
+
+
+ // First of all, here are put some variable declarations
+ // that are useful in the following part:
+ Int_t nparticles; // number of particles
+ // ITS module coordinates [layer = 1, ladder = 2, det = 3] and absolute ID[0] of module [0]
+ TArrayI ID(4);
+ Int_t nmodules, dtype; // Total number of modules and module type (SSD, SPD, SDD)
+ Float_t *x = 0, *y = 0, *z = 0; // Arrays where to store read coords
+ Bool_t *St = 0; // Status of the track (hits only)
+ /*
+ // It's necessary to load the gAlice shared libs
+ // if they aren't already stored in memory...
+ 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;
+ }
+ // }
+
+ // Now is opened the Root input file containing Geometry, Kine and Hits
+ // by default its name must be "galice.root".
+ // When the file is opened, its contens are shown.
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
+ if (!file) file = new TFile(filename);
+ file->ls();
+
+ // Then, the macro gets the AliRun object from file.
+ // If this object is not present, an error occurs
+ // and the execution is stopped.
+ // Anyway, this operation needs some time,
+ // don't worry about an apparent absence of output and actions...
+ cout << "\nSearching in '" << filename << "' for an AliRun object ... " << flush;
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (gAlice)
+ cout << "FOUND!" << endl;
+ else {
+ cout<<"NOT FOUND! The Macro can't continue!" << endl;
+ return 0;
+ }
+ if(!(FileDigits.Data() == filename)){
+ gAlice->SetTreeDFileName(FileDigits);
+ }
+ if(!(FileRec.Data() == filename)){
+ gAlice->SetTreeRFileName(FileRec);
+ }
+
+ // Then, the macro selects the event number specified. Default is 0.
+ nparticles = gAlice->GetEvent(evNum);
+ cout << "\nNumber of particles = " << nparticles << endl;
+ if (!nparticles) {
+ cout << "With no particles I can't do much... Goodbye!" << endl;
+ return 0;
+ }
+
+ // The next step is filling the ITS from the AliRun object.
+ AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
+ ITS->InitModules(-1, nmodules);
+ cout << "Number of ITS modules = " << nmodules << endl;
+ cout << "\nFilling modules (it takes a while, now)..." << flush;
+ ITS->FillModules(0, 0, nmodules, " ", " ");
+ cout << "DONE!" << endl;
+ AliITSgeom *gm = ITS->GetITSgeom();
+ // AliITSDetType *det = ITS->DetType(dtype);
+ // AliITSsegmentation *seg = det->GetSegmentationModel();
+
+ for(;;) {
+
+ // Input phase.
+ // The macro asks if the user wants to put a detector ID[0]
+ // or prefers to input layer, ladder and detector.
+ for (Int_t i = 0; i < 4; i++) ID[i] = -1;
+ Int_t answ;
+ do {
+ cout << "\nSelection modes:" << endl;
+ cout << "1) by layer - ladder - detector numbers" << endl;
+ cout << "2) by unique detector ID" << endl;
+ cout << "0) exit macro" << endl;
+ cout << "\nEnter your choice: ";
+ cin >> answ;
+ } while (answ < 0 || answ > 2);
+ switch (answ) {
+ case 0:
+ // input = 0 ---> EXIT
+ return 0;
+ break;
+ case 1:
+ // input = 1 ---> SELECTION BY COORD
+ do {
+ cout << "\nLayer number [1-6, 0 to exit]: ";
+ cin >> ID[1];
+ if (!ID[1]) return 0;
+ } while (ID[1] < 0 || ID[1] > 6);
+
+ // Detector type: 0 --> SPD, 1 --> SDD, 2 --> SSD.
+ // Layer 1,2 --> 0 / Layer 3,4 --> 1 / Layer 5,6 --> 2
+ dtype = (ID[1] - 1) / 2;
+
+ // Once fixed the layer number, the macro calculates the max number
+ // for ladder and detector from geometry, and accepts only suitable values.
+ do {
+ ID[2] = gm->GetNladders(ID[1]);
+ cout << "Ladder number [1-" << ID[2] << ", 0 to exit]: ";
+ cin >> ID[2];
+ if (!ID[2]) return 0;
+ } while (ID[2] < 0 || ID[2] > gm->GetNladders(ID[1]));
+ do {
+ ID[3] = gm->GetNdetectors(ID[1]);
+ cout << "Detector number [1-" << ID[3] << ", 0 to exit]: ";
+ cin >> ID[3];
+ if (!ID[3]) return 0;
+ } while (ID[3] < 0 || ID[3] > gm->GetNdetectors(ID[1]));
+ break;
+ case 2:
+ // input = 2 ---> SELECTION BY ID[0]
+ do {
+ ID[0] = gm->GetIndexMax();
+ cout << "\n Detector ID number [0-" << ID[0]-1 << ", -1 to exit]: ";
+ cin >> ID[0];
+ ID[0]++;
+ if (ID[0] == -1) return 0;
+ } while (ID[0] < 0 || ID[0] > gm->GetIndexMax());
+ break;
+ };
+
+ if (ID[0] == -1)
+ // If ID[0] = -1 the chioce was by coords, so it's necessary to assign the ID:
+ ID[0] = gm->GetModuleIndex(ID[1], ID[2], ID[3]);
+ else {
+ // Else we must get the layer, ladder and detector by the ID;
+ // ...but first we must remember that the ID goest from 0 to NModules - 1!
+ ID[0]--;
+ ID[1] = ITS->GetModule(ID[0])->GetLayer();
+ ID[2] = ITS->GetModule(ID[0])->GetLadder();
+ ID[3] = ITS->GetModule(ID[0])->GetDet();
+ }
+
+ // Defines the histograms inside the `for' cycle, so they are destroyed at the end
+ // of every read sequqnce, in order to mek another withour segmentation faults
+ Text_t msg[250];
+ Float_t xm = 0.0, zm = 0.0;
+ dtype = (ID[1] - 1) / 2;
+ switch (dtype) {
+ case 0: xm = 1.5; zm = 7.0; break;
+ case 1: xm = 7.5; zm = 8.0; break;
+ case 2: xm = 7.5; zm = 4.5; break;
+ }
+ sprintf(msg, "Module index=%d lay=%d lad=%d det=%d", ID[0], ID[1], ID[2], ID[3]);
+ TH2F *hhits = new TH2F("hhits", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of hits
+ TH2F *hrecs = new TH2F("hrecs", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of recpoints
+ TH2F *hdigits = new TH2F("hdigits", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of digits
+
+ cout << endl;
+
+ // Reads hits...
+ Int_t tmpint = ID[0];
+ Int_t hits = GetModuleHits(ITS, tmpint, x, y, z, St);
+ if (!hits) {
+ cout << "No hits in module!" << endl;
+ continue;
+ }
+ else
+ cout << "Hits scanned..." << endl;
+ for (Int_t i = 0; i < hits; i++) if (!St[i]) hhits->Fill(x[i], z[i]);
+
+ // Reads recpoints...
+ Int_t recs = GetModuleRecPoints(ITS, ID[0], x, z, isfastpoints);
+ if (!recs) {
+ cout << "No recpoints in module!" << endl;
+ continue;
+ }
+ else
+ cout << "Recpoints scanned..." << endl;
+ for (Int_t i = 0; i < recs; i++) hrecs->Fill(x[i], z[i]);
+
+ // Reads digits...
+ Int_t digits = GetModuleDigits(ITS, ID[0], dtype, x, z);
+ if (!digits) {
+ cout << "No digits in module!" << endl;
+ //continue;
+ }
+ else
+ cout << "Digits scanned..." << endl;
+ for (Int_t i = 0; i < digits; i++) hdigits->Fill(x[i], z[i]);
+
+ // Draws read data...
+ // HITS -------> red (2) crosses.
+ // DIGITS -----> green (8) boxes.
+ // REC-POINTS -> blue (4) St. Andrew's crosses.
+
+ TCanvas *viewer = new TCanvas("viewer", "Module viewer canvas", 0, 0, 800, 800);
+ viewer->cd();
+
+ hdigits->SetMarkerStyle(25);
+ hdigits->SetMarkerColor(8);
+ hdigits->SetMarkerSize(2);
+ hdigits->SetStats(kFALSE);
+ hdigits->SetXTitle("Local X (cm)");
+ hdigits->SetYTitle("Local Z (cm)");
+ hdigits->Draw();
+
+ hhits->SetMarkerStyle(5);
+ hhits->SetMarkerColor(2);
+ hhits->SetMarkerSize(3);
+ hhits->SetStats(kFALSE);
+ hhits->Draw("same");
+
+ hrecs->SetMarkerStyle(2);
+ hrecs->SetMarkerColor(4);
+ hrecs->SetMarkerSize(3);
+ hrecs->SetStats(kFALSE);
+ hrecs->Draw("same");
+
+ TLegend *legend = new TLegend(0.7, 0.8, 0.9, 0.9);
+ legend->SetMargin(0.2);
+ legend->AddEntry(hhits, "hits","P");
+ legend->AddEntry(hrecs, "recpoints","P");
+ legend->AddEntry(hdigits, "digits","P");
+ legend->SetTextSizePixels(14);
+ legend->Draw();
+
+ viewer->Update();
+
+
+ Text_t fname[250],ans;
+ cout << "Do you want to save the current canvas on a file (y/n) ? ";
+ cin >> ans;
+ if(ans == 'y' || ans == 'Y') {
+ cout << "Enter filename: ";
+ cin >> fname;
+ TString *control = new TString(fname);
+ Bool_t ok=control->Contains(".C") || control->Contains(".root") || control->Contains(".ps") || control->Contains(".eps") || control->Contains(".gif");
+ if(!ok){
+ cout << "File extension is not recognized. The canvas will be saved as Postscript file";
+ strcat(fname,".ps");
+ }
+ viewer->SaveAs(fname);
+ }
+ }
+
+ cout << "Done. Goodbye" << endl;
+ return 0;
+}
+
+Int_t GetModuleHits(TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St) {
+ // First of all, the macro selects the specified module,
+ // then gets the array of hits in it and their number.
+ AliITS *ITS = (AliITS*) its;
+ AliITSmodule *module = ITS->GetModule(ID);
+ TObjArray *hits_array = module->GetHits();
+ Int_t hits_num = hits_array->GetEntriesFast();
+
+ // Now, if this count returns 0, there's nothing to do,
+ // while, if it doesn't, the first thing to do is dimensioning
+ // the coordinate arrays, and then the loop can start.
+ if (!hits_num)
+ return 0;
+ else {
+ if (X) delete [] X;
+ if (Y) delete [] Y;
+ if (Z) delete [] Z;
+ if (St) delete [] St;
+ X = new Float_t[hits_num];
+ Y = new Float_t[hits_num];
+ Z = new Float_t[hits_num];
+ St = new Bool_t[hits_num];
+ }
+
+ for (Int_t j = 0; j < hits_num; j++) {
+ AliITShit *hit = (AliITShit*) hits_array->At(j);
+ X[j] = hit->GetXL();
+ Y[j] = hit->GetYL();
+ Z[j] = hit->GetZL();
+ St[j] = hit->StatusEntering();
+ }
+ return hits_num;
+}
+
+Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z, Int_t isfastpoints) {
+
+ // First of all, the macro selects the specified module,
+ // then gets the array of recpoints in it and their number.
+ AliITS *ITS = (AliITS*) its;
+ TTree *TR = gAlice->TreeR();
+ ITS->ResetRecPoints();
+
+ TClonesArray *recs_array = ITS->RecPoints();
+
+ TBranch *branch = 0;
+ if(TR && recs_array){
+ if(isfastpoints==1){
+ branch = gAlice->TreeR()->GetBranch("ITSRecPointsF");
+ cout<<"using fast points\n";
+ }
+ else {
+ branch = gAlice->TreeR()->GetBranch("ITSRecPoints");
+ }
+ if(branch)branch->SetAddress(&recs_array);
+ }
+
+ branch->GetEvent(ID);
+
+ Int_t recs_num = recs_array->GetEntries();
+
+
+ // Now, if this count returns 0, there's nothing to do,
+ // while, if it doesn't, the first thing to do is dimensioning
+ // the coordinate and energy loss arrays, and then the loop can start.
+ if (!recs_num)
+ return 0;
+ else {
+ if (X) delete [] X;
+ if (Z) delete [] Z;
+ X = new Float_t[recs_num];
+ Z = new Float_t[recs_num];
+ }
+ for(Int_t j = 0; j < recs_num; j++) {
+ AliITSRecPoint *recp = (AliITSRecPoint*)recs_array->At(j);
+ X[j] = recp->GetX();
+ Z[j] = recp->GetZ();
+ }
+ return recs_num;
+}
+
+Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z) {
+
+ // First of all, the macro selects the specified module,
+ // then gets the array of recpoints in it and their number.
+ AliITSdigit *digit;
+ AliITS *ITS = (AliITS*)its;
+ TTree *TD = gAlice->TreeD();
+ ITS->ResetDigits();
+ TD->GetEvent(ID);
+ TClonesArray *digits_array = ITS->DigitsAddress(dtype);
+ AliITSgeom *gm = ITS->GetITSgeom();
+ AliITSDetType *det = ITS->DetType(dtype);
+ AliITSsegmentation *seg = det->GetSegmentationModel();
+ TArrayI ssdone(5000); // used to match p and n side digits of strips
+ TArrayI pair(5000); // as above
+ Int_t digits_num = digits_array->GetEntries();
+ // Now, if this count returns 0, there's nothing to do,
+ // while, if it doesn't, the first thing to do is dimensioning
+ // the coordinate and energy loss arrays, and then the loop can start.
+
+ if(dtype==2){
+ Int_t layer, ladder, detec;
+ gm->GetModuleId(ID,layer,ladder,detec);
+ seg->SetLayer(layer);
+ }
+
+ if (!digits_num)
+ return 0;
+ else {
+ cout << "Digits to scan: " << digits_num << endl;
+ if (X) delete [] X;
+ if (Z) delete [] Z;
+ X = new Float_t[digits_num];
+ Z = new Float_t[digits_num];
+ }
+
+ // Get the coordinates of the module
+ if (dtype == 2) {
+ for (Int_t j = 0; j < digits_num; j++) {
+ ssdone[j] = 0;
+ pair[j] = 0;
+ }
+ }
+ // AliITSdigit *digit;
+ for (Int_t j = 0; j < digits_num; j++) {
+ digit = (AliITSdigit*)digits_array->UncheckedAt(j);
+ Int_t iz=digit->fCoord1; // cell number z
+ Int_t ix=digit->fCoord2; // cell number x
+ // Get local coordinates of the element (microns)
+ // ******************************* PARTE CORRETTA ***************************************
+ if(dtype < 2) {
+ Float_t xx, zz; // aggiunta
+ seg->DetToLocal(ix, iz, xx, zz);
+ X[j] = xx; // aggiunta
+ Z[j] = zz; // aggiunta
+ }
+ // ******************************* FINE PARTE CORRETTA ***************************************
+ else {
+ // SSD: if iz==0 ---> N side; if iz==1 P side
+ if (ssdone[j] == 0) {
+ ssdone[j]=1;
+ pair[j]=-1;
+ Bool_t pside = (iz == 1);
+ Bool_t impaired = kTRUE;
+ Int_t pstrip = 0;
+ Int_t nstrip = 0;
+ if (pside) pstrip = ix; else nstrip = ix;
+ for (Int_t k = 0; k < digits_num; k++) {
+ if (ssdone[k] == 0 && impaired) {
+ AliITSdigitSSD *sdigit=(AliITSdigitSSD*)digits_array->UncheckedAt(k);
+ if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) {
+ ssdone[k]=2;
+ pair[j]=k;
+ if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2;
+ impaired=kFALSE;
+ }
+ }
+ }
+ if (!impaired) seg->GetPadCxz(pstrip, nstrip, X[j], Z[j]);
+ X[j] /= 10000.0; // convert microns to cm
+ Z[j] /= 10000.0; // convert microns to cm
+ }
+ }
+ }
+ return digits_num;
+}
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+TFile* AccessFile(TString inFile="galice.root", TString acctype="R");
+void writeAR(TFile * fin, TFile *fou);
+void AliITSSD2D(TString inFile, TString outFile);
+
+void AliITSSDigits2DigitsDubna(TString inFile= "galiceS.root",
+ TString outFile = "galiceD.root"){
+ // This macro takes SDigits and produces Digits. No merging is done
+ // and only one galice.root file is used.
+ // Dynamically link some shared libs
+ TStopwatch timer;
+
+ if(gAlice){
+ delete gAlice;
+ gAlice = 0;
+ } // end if gAlice
+ cout << "Creating digits from summable digits for the ITS..." << endl;
+ AliITSSD2D(inFile,outFile);
+ timer.Stop();
+ timer.Print();
+}
+//______________________________________________________________________
+void AliITSSD2D(TString inFile, TString outFile){
+ AliRunDigitizer * manager = new AliRunDigitizer(1,1);
+ char ftmp[50];
+ sprintf(ftmp,"%s",inFile.Data());
+ TFile *file0 = AccessFile(ftmp);
+ AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
+ if (!ITS) {
+ cerr<<"AliITSHits2DigitsDefault.C : AliITS object not found on file"
+ << endl;
+ return 3;
+ } // end if !ITS
+ if(!(ITS->GetITSgeom())){
+ cerr << " AliITSgeom not found. Can't digitize with out it." << endl;
+ return 4;
+ } // end if
+
+ // For old files, must change SPD noise.
+ AliITSresponseSPDdubna *resp0 = new AliITSresponseSPDdubna();
+ if(ITS->DetType(0)->GetResponseModel() !=0){
+ delete ((AliITSresponse*)ITS->DetType(0)->GetResponseModel());
+ ITS->DetType(0)->ResponseModel(0);
+ } // end if
+ ITS->DetType(0)->ResponseModel(resp0);
+ AliITSsegmentationSPD *seg0 = (AliITSsegmentationSPD*)ITS->DetType(0)->
+ GetSegmentationModel();
+ AliITSsimulationSPDdubna *sim0 = new AliITSsimulationSPDdubna(seg0,resp0);
+ if(ITS->DetType(0)->GetSimulationModel() !=0){
+ delete ((AliITSsimulation*)ITS->DetType(0)->GetSimulationModel());
+ ITS->DetType(0)->SimulationModel(0);
+ } // end if
+ ITS->DetType(0)->SimulationModel(sim0);
+ manager->SetInputStream(0,ftmp);
+ if(outFile != "")manager->SetOutputFile(outFile);
+ AliITSDigitizer *dITS = new AliITSDigitizer(manager);
+ manager->Exec("");
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(inFile);
+ TFile *file2 = 0;
+ if(outFile != ""){
+ file2 = new TFile(outFile,"UPDATE");
+ writeAR(file,file2);
+ } // end if outFile!=""
+ delete manager;
+ if(file){
+ file->Write();
+ } // end if file
+ if(file2){
+ file2->Close();
+ delete file2;
+ } // end if file2
+}
+//______________________________________________________________________
+TFile * AccessFile(TString FileName, TString acctype){
+
+ // Function used to open the input file and fetch the AliRun object
+
+ TFile *retfil = 0;
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(FileName);
+ if (file) {file->Close(); delete file; file = 0;}
+ if(acctype.Contains("U")){
+ file = new TFile(FileName,"update");
+ }
+ if(acctype.Contains("N") && !file){
+ file = new TFile(FileName,"recreate");
+ }
+ if(!file) file = new TFile(FileName); // default readonly
+ if (!file->IsOpen()) {
+ cerr<<"Can't open "<<FileName<<" !" << endl;
+ return retfil;
+ }
+
+ // Get AliRun object from file or return if not on file
+ if (gAlice) {delete gAlice; gAlice = 0;}
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (!gAlice) {
+ cerr << "AliRun object not found on file"<< endl;
+ return retfil;
+ }
+ return file;
+}
+//______________________________________________________________________
+void writeAR(TFile * fin, TFile *fou) {
+ TDirectory *current = gDirectory;
+ TTree *Te;
+ TTree *TeNew;
+ AliHeader *alhe = new AliHeader();
+ Te = (TTree*)fin->Get("TE");
+ Te->SetBranchAddress("Header",&alhe);
+ Te->SetBranchStatus("*",1);
+ fou->cd();
+ TeNew = Te->CloneTree();
+ TeNew->Write(0,TObject::kOverwrite);
+ gAlice->Write(0,TObject::kOverwrite);
+ current->cd();
+ delete alhe;
+ cout<<"AliRun object written to file\n";
+}
--- /dev/null
+////////////////////////////////////////////////////////////
+// Test macro for AliITSPid class and tracking version V2 //
+// Rev. 25 July 2002 v3-08-03 + Root 3.02/07 + RedHat 6.2 //
+////////////////////////////////////////////////////////////
+void
+AliITSSavePIDV2(Int_t evNumber1=0,Int_t evNumber2=0) {
+ const char *filename="AliITStracksV2.root";
+
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
+ if (!file) file = new TFile(filename);
+
+TFile *fpid = new TFile("AliITStracksV2Pid.root","recreate");
+
+
+ Float_t factor=0;
+ if(gAlice!=0)factor=gAlice->Field()->SolenoidField();
+ if(factor==0.){
+ cout<<" ================ WARNING ====================================\n";
+ cout<<" The default magnetic field value of 0.4 T will be used\n";
+ cout<<" ==============================================================\n";
+ factor=4.; // Default mag. field = 0.4T
+ }
+ else {
+ cout<<"AliITSSavePIDV2.C: Magnetic field is "<<factor/10.<<" T\n";
+ }
+ AliKalmanTrack::SetConvConst(1000/0.299792458/factor);
+
+
+
+ AliITSPid *pid=new AliITSPid(100);
+//
+// Loop over events
+//
+ for (int nev=0; nev<= evNumber2; nev++) {
+ char tname[30];
+ sprintf(tname,"TreeT_ITS_%d;1",nev);
+ TTree *tracktree=(TTree*)file->Get(tname);
+ TBranch *tbranch=tracktree->GetBranch("tracks");
+
+ Int_t nentr=tracktree->GetEntries();
+ cout<<"Found "<<nentr<<" ITS tracks in event No "<<nev<<endl;
+
+ char tpidname[30];
+ sprintf(tpidname,"TreeT%d",nev);
+ AliITStrackV2Pid pidtmp;
+ TTree itspidTree(tpidname,"Tree with PID");
+ AliITStrackV2Pid *outpid=&pidtmp;
+ itspidTree.Branch("pids","AliITStrackV2Pid",&outpid,32000,1);
+ AliITStrackV2 *iotrack=0;
+ for (Int_t i=0; i<nentr; i++) {
+ AliITStrackV2 *iotrack=new AliITStrackV2;
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
+ iotrack->CookdEdx();
+ Float_t signal=iotrack->GetdEdx();
+ iotrack->PropagateTo(3.,0.0028,65.19);
+ iotrack->PropagateToVertex();
+ Double_t xk,par[5]; iotrack->GetExternalParameters(xk,par);
+ Float_t lam=TMath::ATan(par[3]);
+ Float_t pt_1=TMath::Abs(par[4]);
+ Float_t mom=0.;
+ if( (pt_1*TMath::Cos(lam))!=0. ){ mom=1./(pt_1*TMath::Cos(lam)); }else{mom=0.;};
+ Float_t phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
+ if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+ if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+ signal=signal/35.;
+ pidtmp.fSignal=signal;
+ pidtmp.fMom=mom;
+ pidtmp.fPhi=phi;
+ pidtmp.fLam=lam;
+ pidtmp.fGlab=TMath::Abs(iotrack->GetLabel());
+ Int_t pcode=pid->GetPcode(signal,mom);
+ pidtmp.fPcode=pcode;
+ pidtmp.fWpi=pid->GetWpi();
+ pidtmp.fWk=pid->GetWk();
+ pidtmp.fWp=pid->GetWp();
+ //cout<<" pcode,sig,mom="<<pcode<<" "<<signal<<" "<<mom<<endl;
+ //cout<<" wpi,wka,wp="<<pidtmp.fWpi<<" "<<pidtmp.fWk<<" "<<pidtmp.fWp<<endl;
+ itspidTree.Fill();
+ delete iotrack;
+ }// End for i (tracks)
+ cout<<"n ev="<<nev<<endl;
+ fpid->Write();
+ }// End for nev (events)
+ file->Close();
+ cout<<"File AliITStracksV2Pid.root written"<<endl;
+ delete file;
+ cout<<"end of AliITSSavePIDV2.C "<<endl;
+ //return;
+///////////////////////////////////////////////////////
+}
+
+
+
+
+
+
+
+
--- /dev/null
+///////////////////////////////////////////////////////////
+// Test macro for AliITStracksV1Pid.root file //
+// JINR Dubna Jan 2002 //
+///////////////////////////////////////////////////////////
+void
+AliITSScanPIDV1(Int_t evNumber1=0,Int_t evNumber2=0) {
+ //................. Prepare histogramms ................
+ TH2F *qplot = new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ TH2F *qplotP= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ TH2F *qplotKa= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ TH2F *qplotPi= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ TH2F *qplotE= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ qplotP.SetMarkerStyle(8); qplotP.SetMarkerColor(kBlack); qplotP.SetMarkerSize(.3);
+ qplotKa.SetMarkerStyle(8); qplotKa.SetMarkerColor(kRed); qplotKa.SetMarkerSize(.3);
+ qplotPi.SetMarkerStyle(8); qplotPi.SetMarkerColor(kBlue); qplotPi.SetMarkerSize(.3);
+ qplotE.SetMarkerStyle(8); qplotE.SetMarkerColor(kGreen); qplotE.SetMarkerSize(.3);
+ //......................................................
+ TH1F *signal_mip = new TH1F("signal_mip","Signal (mips) for track",100,0.,15.);
+
+//*****************************************************************************************************************************************
+
+ const char *filename="itstracks.root";
+
+ ///////////////// 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(filename);
+ if (!file) file = new TFile(filename);
+
+//
+// Loop over events
+//
+ char tname[30];
+ for (int nev=evNumber1; nev<= evNumber2; nev++) {
+
+
+ // for (int nev=0; nev<= evNumber2; nev++) {
+
+ sprintf(tname,"TreeT%d",nev);
+ TTree *tracktree=(TTree*)file->Get(tname);
+ TBranch *tbranch=tracktree->GetBranch("ITStracks");
+ cout<<" nev = "<<nev<<"\n";
+ //cout<<" open the file \n";
+
+ Int_t nentr=tracktree->GetEntries();
+
+ TObjArray tarray(nentr);
+ // AliITSIOTrack *iotrack=0;
+ printf("nentr %d\n",nentr);
+
+ for (Int_t i=0; i<nentr; i++) {
+ AliITSIOTrack *iotrack=new AliITSIOTrack;
+ // tarray.AddAt(new AliITSIOTrack,i);
+ // iotrack=(AliITSiotrack*)tarray.UncheckedAt(i);
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
+ tarray.AddLast(iotrack);
+ }
+ //file->Close();
+
+ AliITSIOTrack *iotrack;
+ for (Int_t i=0; i<nentr; i++) {
+ AliITSIOTrack *iotrack=new AliITSIOTrack;
+ iotrack=(AliITSIOTrack*)tarray.UncheckedAt(i);
+ if(!iotrack) continue;
+
+ Float_t Px=iotrack->GetPx();
+ Float_t Py=iotrack->GetPy();
+ Float_t Pz=iotrack->GetPz();
+
+ Float_t momentum = TMath::Sqrt(Px*Px+Py*Py+Pz*Pz);
+ Float_t dEdx = iotrack->GetdEdx();
+ Int_t pcode = 211;//iotrack->GetPDG();
+
+ signal_mip->Fill(dEdx);
+
+ if(pcode == 2212) qplotP.Fill(1000*momentum,dEdx);
+ if(pcode == 321) qplotKa.Fill(1000*momentum,dEdx);
+ if(pcode == 211) qplotPi.Fill(1000*momentum,dEdx);
+ if(pcode == 11) qplotE.Fill(1000*momentum,dEdx);
+
+ delete iotrack;
+ }
+
+ } // event loop
+ file->Close();
+
+//*****************************************************************************************************************************************
+
+ //...................... Draw histogramms .................
+ TCanvas *c1 = new TCanvas("PID_test","Scan PID ",200,10,900,700);
+ c1->Divide(2,1);
+ //.........................................................
+ c1->cd(1); gPad->SetFillColor(33);
+ signal_mip->Draw();
+
+ c1->cd(2); //gPad->SetFillColor(33);
+ qplot->Draw();
+ qplotP.Draw("same"); qplotKa.Draw("same"); qplotPi.Draw("same"); qplotE.Draw("same");
+ AliITSPid *pid =new AliITSPid(1000);
+ c1->Range(0,0,1300,10);
+ gStyle->SetLineColor(kRed);
+ gStyle->SetLineWidth(2);
+ TLine *lj[3],*lk[3];
+ for(Int_t j=0;j<3;j++){
+ Float_t x1,x2,y1,y2,xx1,xx2,yy1,yy2;
+ x1=pid->cut[j+1][0]; x2=pid->cut[j+2][0];
+ y1=y2=pid->cut[j+2][2];
+ lj[j]=new TLine(1000*x1,y1,1000*x2,y2); lj[j]->Draw();
+ if(j==0){yy1=10.;}else{yy1=lj[j-1]->GetY1();}
+ yy2=lj[j]->GetY1();
+ xx1=xx2=x1;
+ lk[j]=new TLine(1000*xx1,yy1,1000*xx2,yy2); lk[j]->Draw();
+ }
+ //Draw pions-kaons cuts.
+ TLine *mj[7],*mk[7];
+ for(Int_t j=0;j<7;j++){
+ Float_t x1,x2,y1,y2,xx1,xx2,yy1,yy2;
+ x1=pid->cut[j+2][0]; x2=pid->cut[j+3][0];
+ y1=y2=pid->cut[j+3][5];
+ mj[j]=new TLine(1000*x1,y1,1000*x2,y2); mj[j]->Draw();
+ if(j==0){yy1=10.;}else{yy1=mj[j-1]->GetY1();}
+ yy2=mj[j]->GetY1();
+ xx1=xx2=x1;
+ mk[j]=new TLine(1000*xx1,yy1,1000*xx2,yy2); mk[j]->Draw();
+ }
+ cout<<"End of file AliITStracksV1Pid.root "<<endl;
+ return;
+}
+
--- /dev/null
+///////////////////////////////////////////////////////////
+// Test macro for AliITStracksV2Pid.root file //
+// JINR Dubna Jan 2002 //
+///////////////////////////////////////////////////////////
+void
+AliITSScanPIDV2(Int_t evNumber1=0,Int_t evNumber2=0) {
+ //................. Prepare histogramms ................
+ TH2F *qplot = new TH2F("Qtrm","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ TH2F *qplotP= new TH2F("QtrmP","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ TH2F *qplotKa= new TH2F("QtrmKa","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ TH2F *qplotPi= new TH2F("QtrmPi","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ TH2F *qplotE= new TH2F("QtrmE","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ qplotP.SetMarkerStyle(8); qplotP.SetMarkerColor(kBlue); qplotP.SetMarkerSize(.3);
+ qplotKa.SetMarkerStyle(8); qplotKa.SetMarkerColor(kRed); qplotKa.SetMarkerSize(.3);
+ qplotPi.SetMarkerStyle(8); qplotPi.SetMarkerColor(kBlack); qplotPi.SetMarkerSize(.3);
+ qplotE.SetMarkerStyle(8); qplotE.SetMarkerColor(kGreen); qplotE.SetMarkerSize(.3);
+ //......................................................
+ TH1F *signal_mip = new TH1F("signal_mip","Signal (mips) for track",100,0.,15.);
+
+TFile *fpid = new TFile("AliITStracksV2Pid.root","read");
+fpid->ls();
+//
+// Loop over events
+//
+for (int nev=0; nev<= evNumber2; nev++) {
+ char tpidname[30];
+ sprintf(tpidname,"TreeT%d",nev);
+ TTree *tracktree=(TTree*)fpid->Get(tpidname);
+ TBranch *tbranch=tracktree->GetBranch("pids");
+
+ Int_t nentr=tracktree->GetEntries();
+ cout<<"Found PID for "<<nentr<<" ITS V2 tracks on "<<tpidname<<endl;
+
+ AliITStrackV2Pid *iopid=0;
+for(Int_t ii=0;ii<nentr;ii++)
+ {
+ AliITStrackV2Pid *iopid=new AliITStrackV2Pid;
+ tbranch->SetAddress(&iopid);
+ tracktree->GetEvent(ii);
+
+ signal_mip->Fill(iopid->fSignal);
+
+ if(iopid->fPcode ==2212)qplotP.Fill(iopid->fMom,iopid->fSignal);
+ if(iopid->fPcode == 321)qplotKa.Fill(iopid->fMom,iopid->fSignal );
+ if(iopid->fPcode == 211)qplotPi.Fill(iopid->fMom,iopid->fSignal );
+ if(iopid->fPcode == 11)qplotE.Fill(iopid->fMom,iopid->fSignal );
+ /*
+
+if( (iopid->fWp<0.10)||(iopid->fWk<0.0)||(iopid->fWpi<0.0) ){
+ cout<<"PID pcode,fsignal,fmom= "<<iopid->fPcode<<","<<iopid->fSignal<<","<<iopid->fMom<<endl;
+ cout<<"wpi,wka,wp="<<iopid->fWpi<<" "<<iopid->fWk<<" "<<iopid->fWp<<endl;
+ }
+ */
+ delete iopid;
+ }// Enf for ii (tracks)
+ }// End for nev (events)
+ fpid->Close();
+ //...................... Draw histogramms .................
+ TCanvas *c1 = new TCanvas("PID_test","Scan PID ",200,10,900,700);
+ c1->Divide(2,1);
+ //.........................................................
+ c1->cd(1); gPad->SetFillColor(33);
+ signal_mip->Draw();
+
+ c1->cd(2); //gPad->SetFillColor(33);
+ qplot->Draw();
+ qplotP.Draw("same"); qplotKa.Draw("same"); qplotPi.Draw("same"); qplotE.Draw("same");
+
+ AliITSPid *pid =new AliITSPid(100);
+ fcutka.Draw("same"); fcutpr.Draw("same");
+ c1->Print("ITSPIDplot.ps");
+
+ cout<<"End of file AliITStracksV2Pid.root "<<endl;
+ return;
+}
+
--- /dev/null
+Int_t AliITSStoreFindableTracks
+(Int_t nMinClusters = 5, const Text_t *evname = "galice", Int_t evnum = 0)
+{
+ gSystem->SetIncludePath("-I- -I$ALICE_ROOT/ITS -I$ALICE_ROOT/STEER");
+
+ gROOT->LoadMacro("$ALICE_ROOT/ITS/AliITSStoreFindableTracksCompiled.C+");
+ AliITSStoreFindableTracksCompiled(nMinClusters, evname, evnum);
+}
--- /dev/null
+#include <fstream.h>
+
+#include <TClassTable.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TString.h>
+#include <TClonesArray.h>
+#include <TParticle.h>
+
+#include "AliRun.h"
+#include "AliITS.h"
+#include "AliITSgeom.h"
+#include "AliITSRecPoint.h"
+
+Int_t AliITSStoreFindableTracksCompiled
+(Int_t nMinClusters = 5, const Text_t *evname = "galice", Int_t evnum = 0)
+{
+ // Make sure that ALICE objects are loaded
+ if (gAlice) {
+ delete gAlice;
+ gAlice = 0;
+ }
+
+ // Define the names of all involved files
+ TString strEventFile(evname);
+ TString strOutputFile(evname);
+ strEventFile.Append(".root");
+ strOutputFile.Append("_tracks_");
+ strOutputFile += evnum;
+ strOutputFile.Append(".root");
+
+ // Connect the Root Galice file containing Geometry, Kine and Hits
+ TFile *fileEvent = (TFile*)gROOT->GetListOfFiles()->FindObject(strEventFile);
+ if (!fileEvent) fileEvent = new TFile(strEventFile,"UPDATE");
+
+ // Get AliRun object from file
+ gAlice = (AliRun*)fileEvent->Get("gAlice");
+ if (gAlice) cout << "OK, found an AliRun object in file" << endl;
+
+ // Get ITS related objects and data
+ AliITS* ITS =(AliITS *)gAlice->GetDetector("ITS");
+ if (!ITS) {
+ cerr << "ITS object not found!" << endl;
+ return 1;
+ }
+ AliITSgeom *geometry = ITS->GetITSgeom();
+ if (!geometry) {
+ cerr << "ITS geometry object not found!" << endl;
+ return 2;
+ }
+
+ // Count the number of modules per layer
+ Int_t nLadders, nDetectors, mod_min[6], mod_max[6];
+ for(Int_t i = 0; i < 6; i++) {
+ nLadders = geometry->GetNladders(i + 1);
+ nDetectors = geometry->GetNdetectors(i + 1);
+ mod_min[i] = geometry->GetModuleIndex(i + 1, 1, 1);
+ mod_max[i] = geometry->GetModuleIndex(i + 1, nLadders, nDetectors);
+ }
+
+ // Load event and ITS recpoints
+ Int_t nParticles = gAlice->GetEvent(evnum);
+ cout << "Event number: " << evnum << endl;
+ cout << "# particles : " << nParticles <<endl;
+ if (nParticles <= 0) {
+ cerr << "Can't have <= 0 particles!" << endl;
+ return 3;
+ }
+ AliITSRecPoint *recp = 0;
+ TClonesArray *recPoints = ITS->RecPoints();
+ TObjArray *particles = gAlice->Particles();
+ Int_t nTracks = gAlice->GetNtrack(); //FCA correction
+ Bool_t *hitITSLayer[6];
+ for (Int_t i = 0; i < 6; i++) {
+ hitITSLayer[i] = new Bool_t[nTracks];
+ for (Int_t j = 0; j < nTracks; j++) hitITSLayer[i][j] = kFALSE;
+ }
+
+ // Load recpoints in event
+ TTree *TR = gAlice->TreeR();
+ if (!TR) {
+ cerr << "TreeR object not found!" << endl;
+ return 4;
+ }
+
+ // Scan recpoints and define findable tracks
+ Int_t nModules = (Int_t)TR->GetEntries(), nPoints = 0, nEmpty = 0;
+ cout << "Found " << nModules;
+ cout << " entries in the TreeR (must be one per module!)" << endl;
+ for (Int_t layer = 1; layer <= 6; layer++) {
+ for (Int_t mod = mod_min[layer - 1]; mod <= mod_max[layer - 1]; mod++) {
+ ITS->ResetRecPoints();
+ TR->GetEntry(mod);
+ nPoints = recPoints->GetEntries();
+ if(!nPoints) {
+ nEmpty++;
+ continue;
+ }
+ for (Int_t point = 0; point < nPoints; point++) {
+ recp = (AliITSRecPoint*)recPoints->UncheckedAt(point);
+ for (Int_t it = 0; it < 3; it++) {
+ Int_t track = recp->GetLabel(it);
+ if(track < 0) continue;
+ if(track > nTracks) {
+ cout << "Found track index " << track;
+ cout << " whilw gAlice->GetNtrack() = " << nTracks << endl;
+ continue;
+ }
+ hitITSLayer[layer - 1][track] = kTRUE;
+ } // loop over recpoint labels
+ } //loop over points
+ } //loop over modules
+ } //loop over layers
+ cout << "Found " << nEmpty << " empty modules" << endl;
+
+ // Scan the file of tracks in TPC to retrieve the findable TPC tracks
+ TString strLabelsTPC;
+ Int_t label, pdg_code, nFindablesTPC = 0;
+ Double_t dummy;
+ ifstream tpc("good_tracks_tpc");
+ while (tpc >> label >> pdg_code) {
+ for (Int_t i = 0; i < 6; i++) tpc >> dummy;
+ nFindablesTPC++;
+ strLabelsTPC.Append(Form("[%d]", label));
+ }
+
+ // Define the TTree with tracks data by means of a set of variables
+ Int_t nFindablesITS = 0, nFindablesITSTPC = 0;
+ Int_t nhits, tpc_ok, mother, entry = 0;
+ Double_t vx, vy, vz;
+ Double_t px, py, pz, pt;
+
+ TTree *tree = new TTree("Tracks", "Findable tracks in ITS");
+
+ tree->Branch("vx", &vx, "vx/D");
+ tree->Branch("vy", &vy, "vy/D");
+ tree->Branch("vz", &vz, "vz/D");
+ tree->Branch("px", &px, "px/D");
+ tree->Branch("py", &py, "py/D");
+ tree->Branch("pz", &pz, "pz/D");
+ tree->Branch("pt", &pt, "pt/D");
+ tree->Branch("label", &label, "label/I");
+ tree->Branch("entry", &entry, "entry/I");
+ tree->Branch("mother", &mother, "mother/I");
+ tree->Branch("pdg_code", &pdg_code, "pdg_code/I");
+ tree->Branch("nhits", &nhits, "nhits/I");
+ tree->Branch("tpc_ok", &tpc_ok, "tpc_ok/I");
+
+ // Fill the tree
+ cout << endl;
+ TParticle *p = 0;
+ for (Int_t i = 0; i < nTracks; i++) {
+ nhits = 0;
+ for (Int_t j = 0; j < 6; j++) if (hitITSLayer[j][i]) nhits++;
+ if (nhits < nMinClusters) continue;
+ p = gAlice->Particle(i);
+ px = p->Px();
+ py = p->Py();
+ pz = p->Pz();
+ pt = p->Pt();
+ vx = p->Vx();
+ vy = p->Vy();
+ vz = p->Vz();
+ mother = p->GetFirstMother();
+ cout << "Track " << i << " stored\r" << flush;
+ tpc_ok = (strLabelsTPC.Contains(Form("[%d]", i)));
+ pdg_code = p->GetPdgCode();
+ label = i;
+ nFindablesITS++;
+ if (tpc_ok) nFindablesITSTPC++;
+ tree->Fill();
+ entry++;
+ }
+
+ // Save into a file
+ TFile *fileOutput = new TFile(strOutputFile, "recreate");
+ tree->Write();
+ fileOutput->Close();
+
+ cout << "# findable tracks in TPC : " << nFindablesTPC << endl;
+ cout << "# findable tracks in ITS : " << nFindablesITS << endl;
+ cout << "# findable tracks in ITS+TPC : " << nFindablesITSTPC << endl;
+
+ return 0;
+}
--- /dev/null
+#ifndef __CINT__
+#include "iostream.h"
+#endif
+
+Bool_t TPCSortTracks(Int_t event = 0)
+{
+ TFile *fileTracks = TFile::Open("AliTPCtracks.root");
+ TFile *fileClusters = TFile::Open("AliTPCclusters.root");
+ TFile *fileEvent = TFile::Open("galice.root");
+
+ // get TPC parameterization
+ AliTPCParam *param=(AliTPCParam *)fileEvent->Get("75x40_100x60_150x60");
+ if (!param) {
+ cerr << "(TPCSortTracks) ERROR: TPC parameters have not been found!" << endl;
+ return kFALSE;
+ }
+
+ // read and sort tracks
+ Int_t i;
+ TSortedList tracks_list;
+ AliTPCtrack *iotrack = 0;
+ TTree *tracktree = (TTree*)fileTracks->Get(Form("TreeT_TPC_%d", event));
+ Int_t nentr = (Int_t)tracktree->GetEntries();
+ for (i = 0; i < nentr; i++) {
+ iotrack = new AliTPCtrack;
+ tracktree->SetBranchAddress("tracks", &iotrack);
+ tracktree->GetEvent(i);
+ tracks_list.Add(iotrack);
+ }
+ delete tracktree;
+
+ // assign to each track its GEANT label
+ fileClusters->cd();
+ AliTPCtracker *tracker = new AliTPCtracker(param, event);
+ tracker->LoadInnerSectors();
+ tracker->LoadOuterSectors();
+ TListIter iter(&tracks_list);
+ for (i = 0; i < nentr; i++) {
+ iotrack = (AliTPCtrack*)iter.Next();
+ if (!iotrack) {
+ cerr << "(TPCSortTracks) WARNING: Track no. " << i << " is NULL!!!" << endl;
+ continue;
+ }
+ tracker->CookLabel(iotrack, 0.1);
+ }
+ delete tracker;
+
+ // create the new TTree of TPC tracks sorted w.r. to Pt
+ tracktree = new TTree(Form("TreeT_TPC_%d", event),"Tree with TPC tracks sorted w.r to pt");
+ tracktree->Branch("tracks", "AliTPCtrack", &iotrack, 32000, 0);
+ iter.Reset();
+ for (i = 0; i < nentr; i++) {
+ iotrack = (AliTPCtrack*)iter.Next();
+ if (!iotrack) {
+ cerr << "(TPCSortTracks) WARNING: Track no. " << i << " is NULL!!!" << endl;
+ continue;
+ }
+ tracktree->Fill();
+ }
+
+ // save the new tree into new file
+ TFile *fileOutput = TFile::Open("AliTPCtracksSorted.root","recreate");
+ tracktree->Write();
+ fileOutput->Close();
+ fileEvent->Close();
+ fileClusters->Close();
+ fileTracks->Close();
+
+ return kTRUE;
+}
+
+
+void AliITSTrackingV1(Int_t evNumber1=0,Int_t evNumber2=0, Int_t min_t=-1, Int_t max_t=0,Bool_t flagvert=1, Bool_t realmass=0, const char *filename="galice.root") {
+
+ ///////////////// 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(filename);
+ if (!file) file = new TFile("galice.root","UPDATE");
+ //if (!file) file = new TFile(filename);
+
+// 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");
+ }
+
+ AliKalmanTrack::SetMagneticField(gAlice->Field()->SolenoidField() / 10.0);
+
+
+ cout << "Sorting TPC tracks w.r. to transverse momentum...";
+ Bool_t success_sorting = TPCSortTracks();
+ if (success_sorting) {
+ cout << "DONE!" << endl;
+ }
+ else {
+ cout << "Some error occurred..." << endl;
+ return 1;
+ }
+
+ AliITS* IITTSS =(AliITS *)gAlice->GetDetector("ITS");
+ if (!IITTSS) return;
+
+//
+// Loop over events
+//
+ Int_t Nh=0;
+ Int_t Nh1=0;
+ for (Int_t nev=0; nev<= evNumber2; nev++) {
+ AliITSTrackerV1 *ITStracker = new AliITSTrackerV1(IITTSS,nev,flagvert);
+ Int_t nparticles = gAlice->GetEvent(nev);
+ cout << "nev " << nev <<endl;
+ cout << "nparticles " << nparticles <<endl;
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+
+ TTree *TR=gAlice->TreeR();
+ Int_t nent=TR->GetEntries();
+ //printf("Found %d entries in the TreeR (must be one per module per event!)\n",nent);
+
+
+ TStopwatch timer;
+
+ timer.Start();
+ ITStracker->DoTracking(nev,min_t,max_t,file,realmass);
+ timer.Stop(); timer.Print();
+ AliITSgeom *g1 = IITTSS->GetITSgeom();
+ Int_t NumOfModules = g1->GetIndexMax();
+ ITStracker->DelMatrix(NumOfModules);
+ delete ITStracker;
+ } // event loop
+ file->Close();
+}
+
--- /dev/null
+////////////////////////////////////////////////////////////////////////
+//
+// AliITSTrackingV2.C
+//
+// date: 18.03.2003
+// author: Thomas Kuhr based on AliTPCTracking.C, AliITSFindClustersV2.C
+// and AliITSFindTracksV2.C
+// version: 1.0
+// description:
+// reconstructs of tracks in ITS in the following steps:
+// ITS cluster finding
+// ITS track finding
+// input parameters:
+// Int_t nEvents ... nr of events to process
+// Int_t firstEventNr ... first event number (starts from 0)
+// char* fileNameHits ... name of file with hits
+// char* fileNameITSDigits .. name of file with ITS digits
+// char* fileNameITSTracks .. name of file with TPC tracks
+// char* fileNameITSClusters .. name of file with ITS clusters (output)
+// char* fileNameITSTracks .. name of file with ITS tracks (output)
+//
+////////////////////////////////////////////////////////////////////////
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "Riostream.h"
+#include "TFile.h"
+#include "TBenchmark.h"
+#include "AliITS.h"
+#include "AliITSgeom.h"
+#include "AliITSclustererV2.h"
+#include "AliITStrackerV2.h"
+#include "AliRun.h"
+#endif
+
+Int_t gDEBUG = 2;
+
+Int_t AliITSTrackingV2(Int_t nEvents=1, Int_t firstEvent=0,
+ const char* fileNameHits="galice.root",
+ const char* fileNameITSDigits="its.digits.root",
+ const char* fileNameTPCTracks="tpc.tracks.root",
+ const char* fileNameITSClusters="its.clusters.root",
+ const char* fileNameITSTracks="its.tracks.root");
+
+Int_t ITSFindClustersV2(const char* fileNameITSDigits,
+ const char* fileNameITSClusters,
+ Int_t n, Int_t first);
+Int_t ITSFindTracksV2(const char* fileNameTPCTracks,
+ const char* fileNameITSClusters,
+ const char* fileNameITSTracks,
+ Int_t n, Int_t first);
+
+////////////////////////////////////////////////////////////////////////
+Int_t AliITSTrackingV2(Int_t nEvents, Int_t firstEvent,
+ const char* fileNameHits,
+ const char* fileNameITSDigits,
+ const char* fileNameTPCTracks,
+ const char* fileNameITSClusters,
+ const char* fileNameITSTracks) {
+
+ AliTracker::SetFieldFactor(fileNameHits,kFALSE);
+
+// find clusters
+
+ if (fileNameITSDigits && fileNameITSClusters) {
+ if(ITSFindClustersV2(fileNameITSDigits,fileNameITSClusters,nEvents,firstEvent)) {
+ cerr << "ITS clustering failed \n";
+ return 1;
+ }
+ }
+
+// find tracks
+
+ if (fileNameTPCTracks && fileNameITSClusters && fileNameITSTracks) {
+ if(ITSFindTracksV2(fileNameTPCTracks,fileNameITSClusters,fileNameITSTracks,nEvents,firstEvent)) {
+ cerr << "ITS tracking failed \n";
+ return 2;
+ }
+ }
+
+ return 0;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t ITSFindClustersV2(const char* fileNameITSDigits, const char* fileNameITSClusters, Int_t nEvents, Int_t firstEvent) {
+//
+// create ITS clusters, store them in the file fileNameITSClusters
+// gAlice object must be in memory
+
+ const char *name="ITSFindClustersV2";
+ if (gDEBUG>1) cout<<name<<" starts...\n";
+ if (gDEBUG>1) gBenchmark->Start(name);
+
+ TFile *in =TFile::Open(fileNameITSDigits);
+ if (!in->IsOpen()) {
+ cerr<<"Can't open file "<<fileNameITSDigits<<endl;
+ return 1;
+ }
+ gAlice->SetTreeDFileName(fileNameITSDigits);
+ TFile *out=TFile::Open(fileNameITSClusters,"recreate");
+ if (!out->IsOpen()) {
+ cerr<<"Can't open file "<<fileNameITSClusters<<endl;
+ return 1;
+ }
+
+ AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
+ if (!ITS) { cerr<<"Can't find the ITS !\n"; return 3; }
+
+ AliITSgeom *geom=ITS->GetITSgeom();
+ out->cd();
+ geom->Write();
+ gROOT->cd();
+
+ AliITSclustererV2 clusterer(geom);
+ for (Int_t iEvent = firstEvent; iEvent<firstEvent+nEvents; iEvent++){
+ cout<<"ITSFindClusters: processing event "<<iEvent<<endl;
+ gAlice->GetEvent(iEvent);
+ in->cd();
+ clusterer.SetEvent(iEvent);
+ clusterer.Digits2Clusters(in,out);
+ }
+
+ out->Close();
+ in->Close();
+
+ if (gDEBUG>1) gBenchmark->Show(name);
+ return 0;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t ITSFindTracksV2(const char* fileNameTPCTracks,
+ const char* fileNameITSClusters,
+ const char* fileNameITSTracks,
+ Int_t nEvents, Int_t first) {
+ Int_t rc=0;
+ const char *name="ITSFindTracksV2";
+ if (gDEBUG>1) cout<<name<<" starts...\n";
+ if (gDEBUG>1) gBenchmark->Start(name);
+
+ TFile *out=TFile::Open(fileNameITSTracks,"recreate");
+ if (!out->IsOpen()) {
+ cerr<<"Can't open file "<<fileNameITSTracks<<endl;
+ return 1;
+ }
+ TFile *in =TFile::Open(fileNameTPCTracks);
+ if (!in->IsOpen()) {
+ cerr<<"Can't open file "<<fileNameTPCTracks<<endl;
+ return 1;
+ }
+ TFile *in2 =TFile::Open(fileNameITSClusters);
+ if (!in2->IsOpen()) {
+ cerr<<"Can't open file "<<fileNameITSClusters<<endl;
+ return 1;
+ }
+
+ AliITSgeom *geom=(AliITSgeom*)in2->Get("AliITSgeom");
+ if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
+
+ AliITStrackerV2 tracker(geom);
+ for (Int_t iEvent=first;iEvent<first+nEvents;iEvent++){
+ cout<<"ITSFindTracks -- event "<<iEvent<<endl;
+ tracker.SetEventNumber(iEvent);
+ rc=tracker.Clusters2Tracks(in,out);
+ }
+
+ in->Close();
+ in2->Close();
+ out->Close();
+
+ if (gDEBUG>1) gBenchmark->Show(name);
+ return rc;
+}
+
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+//-- --- standard headers-------------
+#include <Riostream.h>
+//--------Root headers ---------------
+#include <TSystem.h>
+#include <TFile.h>
+#include <TStopwatch.h>
+#include <TObject.h>
+#include <TTree.h>
+//----- AliRoot headers ---------------
+#include "alles.h"
+#include "AliRun.h"
+#include "AliMagF.h"
+#include "AliKalmanTrack.h"
+#include "AliITSVertex.h"
+#include "AliITSVertexer.h"
+#include "AliITSVertexerTracks.h"
+//-------------------------------------
+#endif
+void AliITSVertexerTracksTest(Int_t evFirst=0,Int_t evLast=0,
+ const Char_t *galiceName="galice.root",
+ const Char_t *trksName="AliITStracksV2.root",
+ const Char_t *vtxName="AliITSVertices.root") {
+ /*******************************************************************
+ * *
+ * Test macro for vertexing in pp using tracks. *
+ * Input file must contain trees with AliITStrackV2 objects. *
+ * Output file can be the same file with the tracks *
+ * or another file. *
+ * If the file galice.root is available, B is taken from there, *
+ * otherwise is can be set here "by hand". *
+ * *
+ * Origin: A.Dainese, Padova andrea.dainese@pd.infn.it *
+ *******************************************************************/
+
+ // Look for field value in galice.root
+ Double_t field = 0.4;
+ if(!gSystem->AccessPathName(galiceName,kFileExists)) {
+ TFile *galice = new TFile(galiceName);
+ gAlice = (AliRun*)galice->Get("gAlice");
+ AliMagF *fiel = (AliMagF*)gAlice->Field();
+ field=(Double_t)fiel->SolenoidField()/10.;
+ AliKalmanTrack::SetConvConst(100/0.299792458/field);
+ printf(" B = %3.1f read from gAlice and set\n",field);
+ delete gAlice;
+ gAlice = 0;
+ galice->Close();
+ delete galice;
+ } else {
+ printf(" File galice.root not found: default 0.4 T being used!\n");
+ }
+
+ // Open input and output files
+ TFile *inFile = TFile::Open(trksName);
+ TFile *outFile = TFile::Open(vtxName,"recreate");
+
+ // Create vertexer
+ AliITSVertexerTracks *vertexer =
+ new AliITSVertexerTracks(inFile,outFile,field);
+ vertexer->SetFirstEvent(evFirst);
+ vertexer->SetLastEvent(evLast);
+ vertexer->SetDebug(0);
+ vertexer->SetUseThrustFrame(0);
+ vertexer->PrintStatus();
+ // Find vertices
+ vertexer->FindVertices();
+
+
+ delete vertexer;
+
+
+
+ inFile->Close();
+ outFile->Close();
+ delete inFile;
+ delete outFile;
+
+ return;
+}
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+//-- --- standard headers-------------
+#include <Riostream.h>
+//--------Root headers ---------------
+#include <TSystem.h>
+#include <TFile.h>
+#include <TStopwatch.h>
+#include <TObject.h>
+#include <TTree.h>
+//----- AliRoot headers ---------------
+#include "alles.h"
+#include "AliRun.h"
+#include "AliMagF.h"
+#include "AliKalmanTrack.h"
+#include "AliITSVertex.h"
+#include "AliITSVertexer.h"
+#include "AliITSVertexerTracks.h"
+#include <AliHeader.h>
+#include <AliGenEventHeader.h>
+//-------------------------------------
+#endif
+void AliITSVertexerTracksTest2(Int_t evFirst=0,Int_t evLast=0,
+ const Char_t *galiceName="galice.root",
+ const Char_t *trksName="AliITStracksV2.root",
+ const Char_t *vtxName="AliITSVertices.root") {
+ /*******************************************************************
+ * *
+ * Test macro for vertexing in pp using tracks. *
+ * Input file must contain trees with AliITStrackV2 objects. *
+ * Output file can be the same file with the tracks *
+ * or another file. *
+ * If the file galice.root is available, B is taken from there, *
+ * otherwise is can be set here "by hand". *
+ * *
+ * Origin: A.Dainese, Padova andrea.dainese@pd.infn.it *
+ *******************************************************************/
+
+ // Look for field value in galice.root
+ Double_t field = 0.4;
+ Int_t kDebug = 0;
+ TFile *galice = 0;
+ if(!gSystem->AccessPathName(galiceName,kFileExists)) {
+ galice = new TFile(galiceName);
+ gAlice = (AliRun*)galice->Get("gAlice");
+ AliMagF *fiel = (AliMagF*)gAlice->Field();
+ field=(Double_t)fiel->SolenoidField()/10.;
+ AliKalmanTrack::SetConvConst(100/0.299792458/field);
+ printf(" B = %3.1f read from gAlice and set\n",field);
+ } else {
+ printf(" File galice.root not found: default 0.4 T being used!\n");
+ }
+
+ // Open input and output files
+ TFile *inFile = TFile::Open(trksName);
+ TFile *outFile = TFile::Open(vtxName,"recreate");
+
+ // Create vertexer
+ AliITSVertexerTracks *vertexer =
+ new AliITSVertexerTracks(inFile,outFile,field,0.,0.);
+ vertexer->SetDebug(0);
+ vertexer->SetUseThrustFrame(0);
+ vertexer->PrintStatus();
+
+ AliITSVertex *vert = 0;
+ // Find vertices
+
+ for(Int_t i=evFirst; i<=evLast; i++){
+ if(i%100==0)cout<<"processing event "<<i<<endl;
+ gAlice->GetEvent(i);
+ // The true Z coord. is fetched for comparison
+ AliHeader *header = gAlice->GetHeader();
+ AliGenEventHeader* genEventHeader = header->GenEventHeader();
+ TArrayF primaryVertex(3);
+ genEventHeader->PrimaryVertex(primaryVertex);
+ vert = vertexer->FindVertexForCurrentEvent(i);
+ if(kDebug>0){
+ // Prints the results
+ cout <<"========================================================\n";
+ cout << "Event number: "<<i<<") Z Vertex:"<<endl;
+ if(vert){
+ cout<<"FOUND: "<<vert->GetZv()<<"; "<<vert->GetZRes()<<endl;
+ cout <<" True Z position "<<primaryVertex[2]<<endl;
+ cout<<", diff= "<<(primaryVertex[2]-vert->GetZv())*10000.<<endl;
+ }
+ else {
+ cout<<"NOT FOUND "<<endl;
+ }
+ }
+ if(vert){
+ Double_t pos[3];
+ for(Int_t kk=0;kk<3;kk++)pos[kk]=(Double_t)primaryVertex[kk];
+ vert->SetTruePos(pos);
+ vertexer->WriteCurrentVertex();
+ }
+ }
+
+ delete vertexer;
+
+
+
+ inFile->Close();
+ outFile->Close();
+ delete inFile;
+ delete outFile;
+ galice->Close();
+ delete galice;
+ return;
+}
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+// Standard includes
+#include <iostream.h>
+// Root includes
+#include <TFile.h>
+#include <TObjArray.h>
+#include <TClonesArray.h>
+#include <TTree.h>
+#include <TBranch.h>
+#include <TH1D.h>
+// AliRoot includes
+#include "STEER/AliRun.h"
+#include "ITS/AliITS.h"
+#include "ITS/AliITSRawCluster.h"
+
+#endif
+
+void AliITSplotSPDClusters(const char* filename="galice_80.root"){
+ //------------------------------------------------------------------
+ // This macro will read the TreeC produced during reconstruction and
+ // plot the number and type of clusters found in the SPD.
+ // root [0] .L AliITSplotSPDClusters.C //this loads the macro in memory
+ // root [1] AliITSplotSPDClusters(); //by default process first event
+ // or
+ // root [1] AliITSplotSPDClusters("galice2.root"); // process all events
+ // from the file
+ // galice2.root.
+ // or
+ // root [0] .x AliITSplotSPDClusters.C("galice2.root") // will all of the
+ // events from the
+ // galice2.root
+ //------------------------------------------------------------------
+
+ // delete the current gAlice object, the one from input file will be used
+ if(gAlice){
+ delete gAlice;
+ gAlice = 0;
+ }else{ // Dynamically link some shared libs
+ if(gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ } // end if
+ } // end if gAlice
+
+ // Connect the Root Galice file containing Geometry, Kine and Hits
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
+ if(!file) file = new TFile(filename);
+ // Get AliRun object from file or create it if not on file
+ if(!gAlice) {
+ gAlice = (AliRun*)file->Get("gAlice");
+ if(gAlice) cout << "AliRun object found on file" << endl;
+ if(!gAlice) gAlice = new AliRun("gAlice","SPD Clulster analysis");
+ } // end if !gAlice
+ Int_t nevents = gAlice->GetEventsPerRun();
+
+ // Get pointers to the ITS Alice detectors.
+ AliITS *ITS = (AliITS*) gAlice->GetDetector("ITS");
+ if(ITS==0){
+ cout << "ITS not found. Exiting." << endl;
+ return;
+ } // end if ITS==0
+ AliITSresponseSPDdubna *resp0 = ITS->DetType(0)->GetResponseModel();
+ Float_t diffCoeff= resp0->DistanceOverVoltage();//Get Value of Diffusion Coefficient parameter d/v.
+ diffCoeff *= resp0->Temperature();
+ diffCoeff = TMath::Sqrt(diffCoeff*8.6173376e-5);// K_Boltzman/e_Coulumb
+ char aDiffCoeff[100]; //Character array for sprintf
+
+ sprintf(aDiffCoeff,"Number of SPD Clusters in x, layer 1, DiffCoeff=%f #sqrt{cm}",
+ diffCoeff);
+ TH1D *hclx1 = new TH1D("hclx1",aDiffCoeff,15,0.5,15.5);
+ sprintf(aDiffCoeff,"Number of SPD Clusters in z, layer 1, DiffCoeff=%f #sqrt{cm}",
+ diffCoeff);
+ TH1D *hclz1 = new TH1D("hclz1",aDiffCoeff,5,0.5,5.5);
+ sprintf(aDiffCoeff,"Number of SPD Clusters in x, layer 2, DiffCoeff=%f #sqrt{cm}",
+ diffCoeff);
+ TH1D *hclx2 = new TH1D("hclx2",aDiffCoeff,15,0.5,15.5);
+ sprintf(aDiffCoeff,"Number of SPD Clusters in z, layer 2, DiffCoeff=%f #sqrt{cm}",
+ diffCoeff);
+ TH1D *hclz2 = new TH1D("hclz2",aDiffCoeff,5,0.5,5.5);
+ // Create Arrays with clusters from: Data, Ba/Sa Model, old version of
+ // Dubna
+ Float_t dataX1[9] = {0.493, 0.493, 0.0273, 0.00617, 0.00112, 0.000257,
+ 0.000257, 0.000257, 0.000257};
+ //clusters from the data in the x-dimension with
+ //standard paramaters
+ Float_t baSaX1[11] = {0.942, 0.0180, 0.0112, 0.00389, 0.00203, 0.000257,
+ 0.001, 0.000257, 0.000257, 0.000257, 0.001};
+ //same for Ba/Sa model
+ Float_t dubnaX1[7] = {0.889, 0.0789, 0.0151, 0.00389, 0.001, 0.001,
+ 0.001};
+ //same for old version of Dubna model
+ Float_t dataX2[9] = {0.482, 0.482, 0.0325, 0.00791, 0.00157, 0.000275,
+ 0.000275, 0.000275, 0.000275};
+ //clusters from the data in the x-dimension with
+ //optimized paramaters
+ Float_t baSaX2[11] = {0.946, 0.0175, 0.0245, 0.00482, 0.001, 0.001,
+ 0.000275, 0.001, 0.000275, 0.000275, 0.001};
+ //same for Ba/Sa model
+ Float_t dubnaX2[8] = {0.638, 0.325, 0.0275, 0.01, 0.00431, 0.000275,
+ 0.001, 0.001};
+ //same for old version of Dubna model
+ Float_t dataZ1[4] = {0.889, 0.0624, 0.000257, 0.000257};
+ //clusters from the data in the z-dimension with
+ //standard paramaters
+ Float_t baSaZ1[3] = {1., 0.0160, 0.000257}; //same for Ba/Sa model
+ Float_t dubnaZ1[3] = {0.889, 0.0180, 0.000257};
+ //same for old version of Dubna model
+ Float_t dataZ2[4] = {0.889, 0.0624, 0.000257, 0.000257};
+ //clusters from the data in the z-dimension with
+ //optimized paramaters
+ Float_t baSaZ2[4] = {1., 0.0160, 0.00215, 0.000257}; //same for Ba/Sa model
+ Float_t dubnaZ2[3] = {0.889, 0.0412, 0.000257};
+ //same for old version of Dubna model
+
+ TH1F *hdataX1 = new TH1F("hdataX1","Number of SPD Clusters in x, layer 1",
+ 9, 0.5, 9.5);
+ hdataX1->SetMarkerColor(2);
+ hdataX1->SetMarkerStyle(20);
+ hdataX1->SetMarkerSize(0.7);
+ TH1F *hbaSaX1 = new TH1F("hbaSaX1", "Ba/Sa", 11, 0.5, 11.5);
+ TH1F *hdubnaX1 = new TH1F("hdubnaX1", "Old Dubna", 7, 0.5, 7.5);
+ TH1F *hdataX2 = new TH1F("hdataX2","Number of SPD Clusters in x, layer 2",
+ 9, 0.5, 9.5);
+ hdataX2->SetMarkerColor(2);
+ hdataX2->SetMarkerStyle(20);
+ hdataX2->SetMarkerSize(0.7);
+ TH1F *hbaSaX2 = new TH1F("hbaSaX2", "Ba/Sa", 11, 0.5, 11.5);
+ TH1F *hdubnaX2 = new TH1F("hdubnaX2", "Old Dubna", 8, 0.5, 8.5);
+ TH1F *hdataZ1 = new TH1F("hdataZ1","Number of SPD Clusters in z, layer 1",
+ 4, 0.5, 4.5);
+ hdataZ1->SetMarkerColor(2);
+ hdataZ1->SetMarkerStyle(20);
+ hdataZ1->SetMarkerSize(0.7);
+ TH1F *hbaSaZ1 = new TH1F("hbaSaZ1", "Ba/Sa", 3, 0.5, 3.5);
+ TH1F *hdubnaZ1 = new TH1F("hdubnaZ1", "Old Dubna", 3, 0.5, 3.5);
+ TH1F *hdataZ2 = new TH1F("hdataZ2","Number of SPD Clusters in z, layer 2",
+ 4, 0.5, 4.5);
+ hdataZ2->SetMarkerColor(2);
+ hdataZ2->SetMarkerStyle(20);
+ hdataZ2->SetMarkerSize(0.7);
+ TH1F *hbaSaZ2 = new TH1F("hbaSaZ2", "Ba/Sa", 4, 0.5, 4.5);
+ TH1F *hdubnaZ2 = new TH1F("hdubnaZ2", "Old Dubna", 3, 0.5, 3.5);
+
+ Int_t j = 0; //j will be a counter for the loops to fill the histograms
+ //with the values for clusters for the Data, the Ba/Sa model
+ //and the old Dubna model
+ for(j=0; j<9; j++){
+ hdataX1->SetBinContent((j+1), (Double_t)dataX1[j]);
+ hdataX1->SetBinError((j+1), 0.00001);
+ }
+ for(j=0; j<11; j++) hbaSaX1->Fill((Float_t)j +0.5, baSaX1[j]);
+ for(j=0; j<7; j++) hdubnaX1->Fill((Float_t)j +0.5, dubnaX1[j]);
+ for(j=0; j<9; j++){
+ hdataX2->SetBinContent((j+1), (Double_t)dataX2[j]);
+ hdataX2->SetBinError((j+1), 0.00001);
+ }
+ for(j=0; j<11; j++) hbaSaX2->Fill((Float_t)j +0.5, baSaX2[j]);
+ for(j=0; j<8; j++) hdubnaX2->Fill((Float_t)j +0.5, dubnaX2[j]);
+ for(j=0; j<4; j++){
+ hdataZ1->SetBinContent((j+1), (Double_t)dataZ1[j]);
+ hdataZ1->SetBinError((j+1), 0.00001);
+ }
+ for(j=0; j<3; j++) hbaSaZ1->Fill((Float_t)j +0.5, baSaZ1[j]);
+ for(j=0; j<3; j++) hdubnaZ1->Fill((Float_t)j +0.5, dubnaZ1[j]);
+ for(j=0; j<4; j++){
+ hdataZ2->SetBinContent((j+1), (Double_t)dataZ2[j]);
+ hdataZ2->SetBinError((j+1), 0.00001);
+ }
+ for(j=0; j<4; j++) hbaSaZ2->Fill((Float_t)j +0.5, baSaZ2[j]);
+ for(j=0; j<3; j++) hdubnaZ2->Fill((Float_t)j +0.5, dubnaZ2[j]);
+
+ TTree *tc = 0;
+ TBranch *spdBranch = 0;
+ TClonesArray *spdClustcA = new TClonesArray("AliITSRawClusterSPD",1000);
+ AliITSRawClusterSPD *spdClust = 0;
+ char tn[20];
+ Int_t evnt,i,n,nc;
+ Float_t nclx = 0,nclz = 0;
+ for(evnt=0;evnt<nevents;evnt++){
+ sprintf(tn,"TreeC%d",evnt);
+ tc = (TTree*) file->Get(tn);
+ spdBranch = tc->GetBranch("ITSClustersSPD");
+ spdBranch->SetAddress(&spdClustcA);
+ n = (Int_t) tc->GetEntries();
+ for(i=0;i<n;i++){
+ tc->GetEvent(i);
+ nc = spdClustcA->GetEntries();
+ if(nc>240) nc = 240;
+ for(j=0;j<nc;j++){
+ spdClust = (AliITSRawClusterSPD*) spdClustcA->At(j);
+ nclx = spdClust->NclX();
+ nclz = spdClust->NclZ();
+ if(spdClust->Module()<80){
+ hclx1->Fill(nclx-0.5);
+ hclz1->Fill(nclz-0.5);
+ }else if(spdClust->Module()<240){
+ hclx2->Fill(nclx-0.5);
+ hclz2->Fill(nclz-0.5);
+ } //endif
+ } // end for j
+ } // end for i
+ spdClustcA->Clear();
+ delete spdBranch; spdBranch = 0;
+ delete tc; tc = 0;
+ } // end for evnt
+
+ //Begin Process of normalizing histograms
+ Double_t integral = (Double_t) hclx1->Integral();
+ if(integral>0.0) hclx1->Scale(1./integral);
+ integral = hclz1->Integral();
+ if(integral>0.0) hclz1->Scale(1./integral);
+ integral = hclx2->Integral();
+ if(integral>0.0) hclx2->Scale(1./integral);
+ integral = hclz2->Integral();
+ if(integral>0.0) hclz2->Scale(1./integral);
+
+ hclx1->SetMinimum(0.000257);
+ hclx1->SetMaximum(1.01);
+ hclz1->SetMinimum(0.000274);
+ hclz1->SetMaximum(1.01);
+ hclx2->SetMinimum(0.000257);
+ hclx2->SetMaximum(1.01);
+ hclz2->SetMinimum(0.000257);
+ hclz2->SetMaximum(1.01);
+
+ sprintf(aDiffCoeff,"SPD Clusters with Diffusion Coefficent=%f",diffCoeff);
+ TCanvas *cSPDclusters = new TCanvas("cSPDclusters",aDiffCoeff,
+ 400,10,600,776);
+ cSPDclusters->Divide(2, 2);
+
+ cSPDclusters->cd(1);
+ cSPDclusters->SetLogy(1);
+ hclx1->Draw("hist");
+ hdataX1->Draw("same e1p");
+ hbaSaX1->SetLineColor(4);
+ hbaSaX1->SetLineStyle(2);
+ hbaSaX1->Draw("same");
+ hdubnaX1->SetLineColor(3);
+ hdubnaX1->SetLineStyle(4);
+ hdubnaX1->Draw("same");
+ TLegend *l1 = new TLegend(0.55,0.65,0.76,0.82);
+ l1->AddEntry(hclx1,"New simulation","l");
+ l1->AddEntry(hdataX1,"Test Beam Results","p");
+ l1->AddEntry(hbaSaX1,"Bari/Selero simulation","l");
+ l1->AddEntry(hdubnaX1,"Dubna simulation","l");
+ l1->Draw();
+
+ cSPDclusters->cd(2);
+ cSPDclusters->SetLogy(1);
+ hclz1->Draw("hist");
+ hdataZ1->Draw("same e1p");
+ hbaSaZ2->SetLineColor(4);
+ hbaSaZ1->SetLineStyle(2);
+ hbaSaZ1->Draw("same");
+ hdubnaZ1->SetLineColor(3);
+ hdubnaZ1->SetLineStyle(4);
+ hdubnaZ1->Draw("same");
+ TLegend *l2 = new TLegend(0.55,0.65,0.76,0.82);
+ l2->AddEntry(hclz1,"New simulation","l");
+ l2->AddEntry(hdataZ1,"Test Beam Results","p");
+ l2->AddEntry(hbaSaZ1,"Bari/Selero simulation","l");
+ l2->AddEntry(hdubnaZ1,"Dubna simulation","l");
+ l2->Draw();
+
+ cSPDclusters->cd(3);
+ cSPDclusters->SetLogy(1);
+ hclx2->Draw("hist");
+ TLegend *l3 = new TLegend(0.55,0.65,0.76,0.82);
+ l3->AddEntry(hclx2,"New simulation","l");
+ l3->Draw();
+/*
+ hdataX2->Draw("e1p");
+ hbaSaX2->SetLineColor(4);
+ hbaSaX2->SetLineStyle(2);
+ hbaSaX2->Draw("same");
+ hdubnaX2->SetLineColor(3);
+ hdubnaX2->SetLineStyle(4);
+ hdubnaX2->Draw("same");
+*/
+ cSPDclusters->cd(4);
+ cSPDclusters->SetLogy(1);
+ hclz2->Draw("hist");
+ TLegend *l4 = new TLegend(0.55,0.65,0.76,0.82);
+ l4->AddEntry(hclz2,"New simulation","l");
+ l4->Draw();
+/*
+ hdataZ2->Draw("e1p");
+ hbaSaZ2->SetLineColor(4);
+ hbaSaZ2->SetLineStyle(2);
+ hbaSaZ2->Draw("same");
+ hdubnaZ2->SetLineColor(3);
+ hdubnaZ2->SetLineStyle(4);
+ hdubnaZ2->Draw("same");
+*/
+ cSPDclusters->Update();
+
+ if(gROOT->IsBatch()){
+ cSPDclusters->Print("SPDClusters.eps");
+ } // end if gROOT->IsBatch()
+
+}
--- /dev/null
+#ifndef __CINT__
+ #include <Riostream.h>
+ #include "AliITSgeom.h"
+ #include "AliITStrackerV2.h"
+
+ #include "TFile.h"
+ #include "TStopwatch.h"
+#endif
+
+Int_t AliITSrefitV2(Int_t nev=1) {
+ cerr<<"Propagating tracks inward through the ITS...\n";
+
+ TFile *in=TFile::Open("AliTPCrefited.root");
+ if (!in->IsOpen()) {cerr<<"Can't open AliTPCrefited.root !\n"; return 1;}
+
+ TFile *out=TFile::Open("AliITStracksV2.root","update");
+ if (!out->IsOpen()) {
+ cerr<<"Can't open AliITStracksV2.root !\n"; return 2;
+ }
+ TFile *file=TFile::Open("AliITSclustersV2.root");
+ if (!file->IsOpen()) {
+ cerr<<"Can't open AliITSclustersV2.root !\n";return 3;
+ }
+ AliITSgeom *geom=(AliITSgeom*)file->Get("AliITSgeom");
+
+ Int_t rc=0;
+ TStopwatch timer;
+ AliITStrackerV2 tracker(geom);
+ for (Int_t i=0; i<nev; i++) {
+ cerr<<"Processing event number : "<<i<<endl;
+ tracker.SetEventNumber(i);
+ rc=tracker.RefitInward(in,out);
+ }
+ timer.Stop(); timer.Print();
+
+ delete geom;
+
+ file->Close();
+ in->Close();
+ out->Close();
+
+ return rc;
+}
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "iostream.h"
+#include "TDatime.h"
+#include "TFile.h"
+#include "TString.h"
+#include "../STEER/AliRun.h"
+#include "../STEER/AliRunDigitizer.h"
+#include "ITS/AliITSDigitizer.h"
+#include "ITS/AliITS.h"
+#include "ITS/AliITSDetType.h"
+#include "ITS/AliITSresponseSDD.h"
+#include "TStopwatch.h"
+
+Bool_t GaliceITSok();
+TFile* AccessFile(TString inFile="galice.root", TString acctype="R");
+void writeAR(TFile * fin, TFile *fou);
+Int_t ChangeITSDefaults(TFile *hitfile,AliITS *ITS,TString opt="");
+
+#endif
+
+//#define DEBUG
+
+Int_t AliITSsd2d(TString df="galice.root",TString sf1="galice.root",
+ TString sf2="",TString opt=""){
+ // Produce ITS Digits from SDigits, with psible merging.
+
+ // Dynamically link some shared libs
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ } // end if
+ gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSstandard.C");
+
+ TFile *sfp1 = 0, *sfp2 = 0, *dfp = 0;
+ if(!GaliceITSok()){
+ // gAlice not define. Must open a file and read it in.
+ if(df.CompareTo(sf1) == 0 || df.CompareTo("") == 0) {
+ // Input file = output file
+ sfp1 = AccessFile(sf1,"U"); // input file open for update.
+ }else{ // Input file different from output file.
+ sfp1 = AccessFile(sf1,"R"); // input file open as read only
+ // open output file and create TreeR on it
+ dfp = gAlice->InitTreeFile("S",df);
+ } // end if
+ } // end if !GALICEITSOK()
+ AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
+
+ ChangeITSDefaults(sfp1,ITS,opt);
+ // write the AliRun object to the output file if different from input file.
+ if(dfp) {
+ writeAR(sfp1,dfp);
+ dfp->Close(); // Manager will open in update mode.
+ } // end if
+
+ AliRunDigitizer *manager;
+ if(sf2.CompareTo("")==0) { // do not merge.
+ manager = new AliRunDigitizer(1,1);
+ manager->SetInputStream(0,sf1);
+ }else{
+ manager = new AliRunDigitizer(2,1);
+ manager->SetInputStream(0,sf1.Data());
+ manager->SetInputStream(1,sf2.Data());
+ } // end if
+ if (df.CompareTo(sf1) !=0) {
+ manager->SetOutputFile(df);
+ } // end if
+ AliITSDigitizer *dITS = new AliITSDigitizer(manager);
+ if(opt.Contains("ROI")==0) dITS->SetByRegionOfInterestFlag(1);
+
+ TStopwatch timer;
+ timer.Start();
+ manager->Exec("all");
+ timer.Stop();
+ timer.Print();
+ delete manager;
+
+ if(dfp!=0){
+ cout << df << " size =" << df->GetSize() << endl;
+ }else if(sfp1!=0){
+ cout << sf1 << " size =" << sf1->GetSize() << endl;
+ } // end if sf1!=0
+ return 0;
+}
--- /dev/null
+#
+set echo on
+#
+aliroot -b -q grun.C\(2,\"$ALICE_ROOT/ITS/Config_muon.C\"\) >& aliroot_muon.log
+#
+mv galice.root galice_muon.root
+#
+aliroot -q -b $ALICE_ROOT/ITS/AliITSHits2SDigits.C\(\"galice_muon.root\"\,\"galice_muonS.root\"\) >& aliroot_muonS.log
+#
+aliroot -b -q grun.C\(2,\"$ALICE_ROOT/ITS/Config_bg.C\"\) >& aliroot_bg.log
+#
+cp galice.root galice_all.root
+mv galice.root galice_bg.root
+#
+aliroot -q -b $ALICE_ROOT/ITS/AliITSHits2SDigits.C\(\"galice_bg.root\"\,\"galice_bgS.root\"\) >& aliroot_bgS.log
+#
+aliroot -q -b $ALICE_ROOT/ITS/AliITSHits2Digits.C\(\"galice_bg.root\"\,\"galice_D.root\"\) >& aliroot_D.log
+#
+aliroot -q -b $ALICE_ROOT/ITS/AliITSHits2Digits.C\(\"galice_all.root\"\,\"galice_all.root\"\) >& aliroot_allD.log
+#
+aliroot -q -b $ALICE_ROOT/ITS/AliITSMerge.C\(\"galice_mD.root\",\"galice_muonS.root\",\"galice_bgS.root\"\) >& aliroot_m.log
+#
+aliroot -q -b $ALICE_ROOT/ITS/AliITSDigits2RecPoints.C\(\"galice_mD.root\",\"galice_mR.root\"\) >& aliroot_mRec.log
+#
+aliroot -q -b $ALICE_ROOT/ITS/AliITSDigits2RecPoints.C\(\"galice_D.root\",\"galice_R.root\"\) >& aliroot_Rec.log
+#
+aliroot -q -b $ALICE_ROOT/ITS/AliITSDigits2RecPoints.C\(\"galice_all.root\",\"galice_all.root\"\) >& aliroot_allRec.log
+#
+
--- /dev/null
+#!/bin/sh
+# 12 May 2003,Dubna
+# This script processes the following steps for n events (AliRoot v3-09-08,
+# root v3.05/04)
+# (use the parameter n):
+# - the TPC and ITS slow simulation (hits+digits+clusters),
+# - the TPC+ITS tracking,
+# - the TPC and ITS PID
+# - the AliITStracksPidV2.root is created with the following information:
+# fGlab - track number
+# fPcode - particle code after the TPC PID
+# fMom - particle momentum (from the AliTPCtracks.root)
+# fLam - particle lambder angle (from the AliTPCtracks.root)
+# fPhi - particle phi angle (from the AliTPCtracks.root)
+# fSignal- TPC signal (mips)
+# fWpi - the PID weight for the identified pion
+# fWk - the PID weight for the identified kaon
+# fWp - the PID weight for the identified proton
+# (the weghts are tuned for the HIJING events)
+# and the AliITSScanPIDV2.C is the example of macros to read the PID
+# information.
+# See ITSPIDplot.ps after run of this script.
+
+if [ $# = 0 ]; then nev=1; else nev=$1; fi
+if [ $nev = 0 ]; then nev=1; fi
+#
+# delete eventual old files from the last run
+echo "Start simulation for " $nev " event(s)"
+$ALICE_ROOT/ITS/AliITSDeleteOldFiles.sh
+#
+# run the hit generation
+aliroot -q -b "$ALICE_ROOT/macros/grun.C($nev)"
+# digitize TPC
+aliroot -q -b "$ALICE_ROOT/TPC/AliTPCHits2Digits.C($nev)"
+# TPC tracking
+ln -s galice.root digits.root
+aliroot -q -b "$ALICE_ROOT/TPC/AliTPCFindClustersMI.C($nev)"
+aliroot -b <<EOI
+.L $ALICE_ROOT/TPC/AliTPCFindTracksMI.C
+AliTPCFindTracks($nev);
+.q
+EOI
+#
+# digitize ITS
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSHits2SDigits.C"
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSSDigits2Digits.C"
+# ITS-TPC tracking
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSFindClustersV2.C('s',$nev)"
+aliroot -b <<EOI
+TFile *inkin = TFile::Open("galice.root");
+if (!inkin->IsOpen()) cerr<<"Can't open galice.root !\n";
+if (gAlice) {delete gAlice; gAlice=0;}
+
+gAlice = (AliRun*)inkin->Get("gAlice");
+cout<<"AliRun object found on file "<<gAlice<<endl;
+cout<<"!!!! field ="<<gAlice->Field()->SolenoidField()<<endl;
+AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
+inkin->Close();
+.x $ALICE_ROOT/ITS/AliITSFindTracksV2.C($nev);
+.q
+EOI
+#
+# Do the PID procedure for the ITS. The output file AliITStrackV2Pid.root is
+# created with the information for each track:
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSSavePIDV2.C(0,$[$nev-1])"
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSScanPIDV2.C(0,$[$nev-1])"
+
+
+
--- /dev/null
+#!/bin/sh
+
+# delete eventual old files from the last run
+./AliITSDeleteOldFiles.sh
+
+# do everything for the TPC
+aliroot -q -b "$ALICE_ROOT/TPC/AliTPCtest.C"
+
+# create summable digits for the ITS
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSHits2SDigits.C"
+
+# go from summable digits to digits for the ITS
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSSDigits2Digits.C"
+
+# create reconstructed points for the ITS
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSDigits2RecPoints.C"
+
+# do the tracking V1
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSTrackingV1.C"
+
+# prepare results of tracking V1 for comparison
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSStoreFindableTracks.C"
+
+# do ITS tracking V1 comparison
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSComparisonV1.C"
+
+#
+# after all of the above you can run AliITSPlotTracksV1.C macro under
+# aliroot to see plots of the efficiency and track parameter resolution
+# or AliITSDisplayTracksV1.C to display found tracks on top of the
+# ITS detailed geometry
--- /dev/null
+void Config(){
+ // Set Random Number seed
+ // gRandom->SetSeed(12345);
+ // libraries required by geant321
+ gSystem->Load("libgeant321");
+ new TGeant3("C++ Interface to Geant3");
+ if (!gSystem->Getenv("CONFIG_FILE")){
+ cout<<"Config.C: Creating Run Loader ..."<<endl;
+ AliRunLoader *rl = AliRunLoader::Open("galice.root",
+ AliConfig::fgkDefaultEventFolderName,"recreate");
+ if (rl == 0x0){
+ gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
+ return;
+ } // end if rl==0x0
+ rl->SetCompressionLevel(2);
+ rl->SetNumberOfEventsPerFile(1000);
+ gAlice->SetRunLoader(rl);
+ } // end if !gSystem
+ TGeant3 *geant3 = (TGeant3 *) gMC;
+ // Set External decayer
+ AliDecayer *decayer = new AliDecayerPythia();
+ decayer->SetForceDecay(kAll);
+ decayer->Init();
+ gMC->SetExternalDecayer(decayer);
+ //=======================================================================
+ // ******* GEANT STEERING parameters FOR ALICE SIMULATION *******
+ geant3->SetTRIG(1); //Number of events to be processed
+ geant3->SetSWIT(4, 10);
+ geant3->SetDEBU(0, 0, 1);
+ //geant3->SetSWIT(2,2);
+ geant3->SetDCAY(1);
+ geant3->SetPAIR(1);
+ geant3->SetCOMP(1);
+ geant3->SetPHOT(1);
+ geant3->SetPFIS(0);
+ geant3->SetDRAY(0);
+ geant3->SetANNI(1);
+ geant3->SetBREM(1);
+ geant3->SetMUNU(1);
+ geant3->SetCKOV(1);
+ geant3->SetHADR(1);//Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3)
+ geant3->SetLOSS(2);
+ geant3->SetMULS(1);
+ geant3->SetRAYL(1);
+ geant3->SetAUTO(1);//Select automatic STMIN etc... calc. (AUTO 1) or manual (AUTO 0)
+ geant3->SetABAN(0);//Restore 3.16 behaviour for abandoned tracks
+ geant3->SetOPTI(2);//Select optimisation level for GEANT geometry searches (0,1,2)
+ geant3->SetERAN(5.e-7);
+ Float_t cut = 1.e-3; // 1MeV cut by default
+ Float_t tofmax = 1.e10;
+ // GAM ELEC NHAD CHAD MUON EBREM MUHAB EDEL MUDEL MUPA TOFMAX
+ geant3->SetCUTS(cut, cut, cut, cut, cut, cut, cut, cut, cut, cut,
+ tofmax);
+ //=======================================================================
+ // ************* STEERING parameters FOR ALICE SIMULATION **************
+ // --- Specify event type to be tracked through the ALICE setup
+ // --- All positions are in cm, angles in degrees, and P and E in GeV
+ if (gSystem->Getenv("CONFIG_NPARTICLES")){
+ int nParticles = atoi(gSystem->Getenv("CONFIG_NPARTICLES"));
+ }else{
+ int nParticles = 1;
+ } // end if
+ //*********************************************
+ // Example for Moving Particle Gun *
+ //*********************************************
+ AliGenBox *gener = new AliGenBox(nParticles);
+ gener->SetMomentumRange(100.,300.);
+ gener->SetPhiRange(0,0.1);
+ gener->SetThetaRange(0.0, .1);
+ gener->SetOrigin(0.,0.,-50.);
+ //vertex position
+ gener->SetSigma(0.1,0.1,0.0); //Sigma in (X,Y,Z) (cm) on IP position
+ gener->SetPart(211); //GEANT particle type
+ gener->Init();
+ // Activate this line if you want the vertex smearing to happen
+ // track by track
+ //
+ //gener->SetVertexSmear(perTrack);
+ // Field (L3 0.4 T)
+ rootfile->cd();
+ //gAlice->SetField(field);
+
+ Int_t iHALL = 0;
+ Int_t iITS = 1;
+ rl->CdGAFile();
+ //=================== Alice BODY parameters =============================
+ AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
+
+ if (iHALL){
+ //=================== HALL parameters ============================
+ AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
+ } // end if
+ if(iITS) {
+ //=================== ITS parameters ============================
+ AliITSvSPD02 *ITS = new AliITSvSPD02("SPD test beam 2002");
+ }
+ return;
+}
+
+Float_t EtaToTheta(Float_t arg){
+ return (180./TMath::Pi())*2.*atan(exp(-arg));
+}
--- /dev/null
+void Config(){
+ // 7-DEC-2000 09:00
+ // Switch on Transition Radiation simulation. 6/12/00 18:00
+ // iZDC=1 7/12/00 09:00
+ // ThetaRange is (0., 180.). It was (0.28,179.72) 7/12/00 09:00
+ // Theta range given through pseudorapidity limits 22/6/2001
+
+ // Set Random Number seed
+ // gRandom->SetSeed(12345);
+
+ // libraries required by geant321
+ gSystem->Load("libgeant321");
+
+ new TGeant3("C++ Interface to Geant3");
+
+ if (!gSystem->Getenv("CONFIG_FILE"))
+ {
+ TFile *rootfile = new TFile("galice.root", "recreate");
+
+ rootfile->SetCompressionLevel(2);
+ }
+
+ TGeant3 *geant3 = (TGeant3 *) gMC;
+
+ //
+ // Set External decayer
+ AliDecayer *decayer = new AliDecayerPythia();
+
+ decayer->SetForceDecay(kAll);
+ decayer->Init();
+ gMC->SetExternalDecayer(decayer);
+ //
+ //
+ //=======================================================================
+ // ******* GEANT STEERING parameters FOR ALICE SIMULATION *******
+ geant3->SetTRIG(1); //Number of events to be processed
+ geant3->SetSWIT(4, 10);
+ geant3->SetDEBU(0, 0, 1);
+ //geant3->SetSWIT(2,2);
+ geant3->SetDCAY(1);
+ geant3->SetPAIR(1);
+ geant3->SetCOMP(1);
+ geant3->SetPHOT(1);
+ geant3->SetPFIS(0);
+ geant3->SetDRAY(0);
+ geant3->SetANNI(1);
+ geant3->SetBREM(1);
+ geant3->SetMUNU(1);
+ geant3->SetCKOV(1);
+ geant3->SetHADR(1); //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3)
+ geant3->SetLOSS(2);
+ geant3->SetMULS(1);
+ geant3->SetRAYL(1);
+ geant3->SetAUTO(1); //Select automatic STMIN etc... calc. (AUTO 1) or manual (AUTO 0)
+ geant3->SetABAN(0); //Restore 3.16 behaviour for abandoned tracks
+ geant3->SetOPTI(2); //Select optimisation level for GEANT geometry searches (0,1,2)
+ geant3->SetERAN(5.e-7);
+
+ Float_t cut = 1.e-3; // 1MeV cut by default
+ Float_t tofmax = 1.e10;
+
+ // GAM ELEC NHAD CHAD MUON EBREM MUHAB EDEL MUDEL MUPA TOFMAX
+ geant3->SetCUTS(cut, cut, cut, cut, cut, cut, cut, cut, cut, cut,
+ tofmax);
+ //
+ //=======================================================================
+ // ************* STEERING parameters FOR ALICE SIMULATION **************
+ // --- Specify event type to be tracked through the ALICE setup
+ // --- All positions are in cm, angles in degrees, and P and E in GeV
+ if (gSystem->Getenv("CONFIG_NPARTICLES")){
+ int nParticles = atoi(gSystem->Getenv("CONFIG_NPARTICLES"));
+ } else{
+ int nParticles = 500;
+ } // end if
+ AliGenHIJINGpara *gener = new AliGenHIJINGpara(nParticles);
+
+ gener->SetMomentumRange(0, 999);
+ gener->SetPhiRange(0, 360);
+ // Set pseudorapidity range from -8 to 8.
+ Float_t thmin = EtaToTheta(8); // theta min. <---> eta max
+ Float_t thmax = EtaToTheta(-8); // theta max. <---> eta min
+ gener->SetThetaRange(thmin,thmax);
+ gener->SetOrigin(0, 0, 0); //vertex position
+ gener->SetSigma(0, 0, 0); //Sigma in (X,Y,Z) (cm) on IP position
+ gener->Init();
+ //
+ // Activate this line if you want the vertex smearing to happen
+ // track by track
+ //
+ //gener->SetVertexSmear(perTrack);
+
+ gAlice->SetField(-999, 2); //Specify maximum magnetic field in Tesla (neg. ==> default field)
+
+ Int_t iABSO = 1;
+ Int_t iDIPO = 1;
+ Int_t iFMD = 1;
+ Int_t iFRAME = 1;
+ Int_t iHALL = 1;
+ Int_t iITS = 1;
+ Int_t iMAG = 1;
+ Int_t iMUON = 1;
+ Int_t iPHOS = 1;
+ Int_t iPIPE = 1;
+ Int_t iPMD = 1;
+ Int_t iRICH = 1;
+ Int_t iSHIL = 1;
+ Int_t iSTART = 1;
+ Int_t iTOF = 1;
+ Int_t iTPC = 1;
+ Int_t iTRD = 1;
+ Int_t iZDC = 1;
+ Int_t iEMCAL = 1;
+
+ //=================== Alice BODY parameters =============================
+ AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
+
+ if (iMAG){
+ //=================== MAG parameters ============================
+ // --- Start with Magnet since detector layouts may be depending ---
+ // --- on the selected Magnet dimensions ---
+ AliMAG *MAG = new AliMAG("MAG", "Magnet");
+ }
+ if (iABSO){
+ //=================== ABSO parameters ============================
+ AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber");
+ }
+ if (iDIPO){
+ //=================== DIPO parameters ============================
+ AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
+ }
+ if (iHALL){
+ //=================== HALL parameters ============================
+ AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
+ }
+ if (iFRAME){
+ //=================== FRAME parameters ============================
+ AliFRAME *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
+ }
+ if (iSHIL){
+ //=================== SHIL parameters ============================
+ AliSHIL *SHIL = new AliSHILvF("SHIL", "Shielding");
+ }
+ if (iPIPE){
+ //=================== PIPE parameters ============================
+ AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
+ }
+ if(iITS) {
+ //=================== ITS parameters ============================
+ //
+ // As the innermost detector in ALICE, the Inner Tracking System
+ // "impacts" on almost all other detectors. This involves the fact
+ // that the ITS geometry still has several options to be followed
+ //in parallel in order to determine the best set-up which minimizes
+ // the induced background. All the geometries available to date are
+ // described in the following. Read carefully the comments and use
+ // the default version (the only one uncommented) unless you are making
+ // comparisons and you know what you are doing. In this case just
+ // uncomment the ITS geometry you want to use and run Aliroot.
+ // Detailed geometries:
+ AliITSvPPRasymm *ITS = new AliITSvPPRasymm("ITS","New ITS PPR detailed version with asymmetric services");
+ ITS->SetMinorVersion(2); // don't touch this parameter if you're not an ITS developer
+ ITS->SetReadDet(kFALSE); // don't touch this parameter if you're not an ITS developer
+ ITS->SetThicknessDet1(200.); // detector thickness on layer 1 must be in the range [100,300]
+ ITS->SetThicknessDet2(200.); // detector thickness on layer 2 must be in the range [100,300]
+ ITS->SetThicknessChip1(200.); // chip thickness on layer 1 must be in the range [150,300]
+ ITS->SetThicknessChip2(200.);// chip thickness on layer 2 must be in the range [150,300]
+ ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ ITS->SetCoolingFluid(1); // 1 --> water ; 0 --> freon
+ // Geant3 <-> EUCLID conversion
+ // ============================
+ // SetEUCLID is a flag to output (=1) or not to output (=0) both
+ // geometry and media to two ASCII files (called by default
+ // ITSgeometry.euc and ITSgeometry.tme) in a format understandable
+ // to the CAD system EUCLID. The default (=0) means that you dont
+ // want to use this facility.
+ ITS->SetEUCLID(0);
+ }
+ if (iTPC){
+ //============================ TPC parameters =========================
+ // This allows the user to specify sectors for the SLOW (TPC geometry
+ // 2) Simulator. SecAL (SecAU) <0 means that ALL lower (upper)
+ // sectors are specified, any value other than that requires at least
+ // one sector (lower or upper)to be specified!
+ // Reminder: sectors 1-24 are lower sectors ( 1-12->z>0,13-24->z<0)
+ // sectors 25-72 are the upper ones (25-48->z>0,49-72->z<0)
+ // SecLows - number of lower sectors specified (up to 6)
+ // SecUps - number of upper sectors specified (up to 12)
+ // Sens - sensitive strips for the Slow Simulator !!!
+ // This does NOT work if all S or L-sectors are specified, i.e.
+ // if SecAL or SecAU < 0
+ //---------------------------------------------------------------------
+ // gROOT->LoadMacro("SetTPCParam.C");
+ // AliTPCParam *param = SetTPCParam();
+ AliTPC *TPC = new AliTPCv2("TPC", "Default");
+ // All sectors included
+ TPC->SetSecAL(-1);
+ TPC->SetSecAU(-1);
+ }
+ if (iTOF){
+ //=================== TOF parameters ============================
+ AliTOF *TOF = new AliTOFv2("TOF", "normal TOF");
+ }
+ if (iRICH){
+ //=================== RICH parameters ===========================
+ AliRICH *RICH = new AliRICHv3("RICH", "normal RICH");
+ }
+ if (iZDC){
+ //=================== ZDC parameters ============================
+ AliZDC *ZDC = new AliZDCv2("ZDC", "normal ZDC");
+ }
+ if (iTRD){
+ //=================== TRD parameters ============================
+ AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
+ // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 90% Xe +
+ // 10% CO2)
+ TRD->SetGasMix(1);
+ // With hole in front of PHOS
+ TRD->SetPHOShole();
+ // With hole in front of RICH
+ TRD->SetRICHhole();
+ // Switch on TR
+ AliTRDsim *TRDsim = TRD->CreateTR();
+ }
+ if (iFMD){
+ //=================== FMD parameters ============================
+ AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
+ FMD->SetRingsSi1(256);
+ FMD->SetRingsSi2(64);
+ FMD->SetSectorsSi1(20);
+ FMD->SetSectorsSi2(24);
+ }
+ if (iMUON){
+ //=================== MUON parameters ===========================
+ AliMUON *MUON = new AliMUONv1("MUON", "default");
+ }
+ if (iPHOS){
+ //=================== PHOS parameters ===========================
+ AliPHOS *PHOS = new AliPHOSv1("PHOS", "GPS2");
+ }
+ if (iPMD){
+ //=================== PMD parameters ============================
+ AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
+ PMD->SetPAR(1., 1., 0.8, 0.02);
+ PMD->SetIN(6., 18., -580., 27., 27.);
+ PMD->SetGEO(0.0, 0.2, 4.);
+ PMD->SetPadSize(0.8, 1.0, 1.0, 1.5);
+ }
+ if (iEMCAL){
+ //=================== EMCAL parameters ============================
+ AliEMCAL *EMCAL = new AliEMCALv1("EMCAL", "G56_2_55_19_104_14");
+ }
+ if (iSTART){
+ //=================== START parameters ============================
+ AliSTART *START = new AliSTARTv1("START", "START Detector");
+ }
+}
+//----------------------------------------------------------------------
+Float_t EtaToTheta(Float_t arg){
+ return (180./TMath::Pi())*2.*atan(exp(-arg));
+}
--- /dev/null
+void Config(){
+ // 7-DEC-2000 09:00
+ // Switch on Transition Radiation simulation. 6/12/00 18:00
+ // iZDC=1 7/12/00 09:00
+ // ThetaRange is (0., 180.). It was (0.28,179.72) 7/12/00 09:00
+ // Theta range given through pseudorapidity limits 22/6/2001
+
+ // Set Random Number seed
+ // gRandom->SetSeed(12345);
+
+ // libraries required by geant321
+ gSystem->Load("libgeant321");
+
+ new TGeant3("C++ Interface to Geant3");
+
+ if (!gSystem->Getenv("CONFIG_FILE"))
+ {
+ TFile *rootfile = new TFile("galice.root", "recreate");
+
+ rootfile->SetCompressionLevel(2);
+ }
+
+ TGeant3 *geant3 = (TGeant3 *) gMC;
+
+ //
+ // Set External decayer
+ AliDecayer *decayer = new AliDecayerPythia();
+
+ decayer->SetForceDecay(kAll);
+ decayer->Init();
+ gMC->SetExternalDecayer(decayer);
+ //
+ //
+ //=======================================================================
+ // ******* GEANT STEERING parameters FOR ALICE SIMULATION *******
+ geant3->SetTRIG(1); //Number of events to be processed
+ geant3->SetSWIT(4, 10);
+ geant3->SetDEBU(0, 0, 1);
+ //geant3->SetSWIT(2,2);
+ geant3->SetDCAY(1);
+ geant3->SetPAIR(1);
+ geant3->SetCOMP(1);
+ geant3->SetPHOT(1);
+ geant3->SetPFIS(0);
+ geant3->SetDRAY(0);
+ geant3->SetANNI(1);
+ geant3->SetBREM(1);
+ geant3->SetMUNU(1);
+ geant3->SetCKOV(1);
+ geant3->SetHADR(1); //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3)
+ geant3->SetLOSS(2);
+ geant3->SetMULS(1);
+ geant3->SetRAYL(1);
+ geant3->SetAUTO(1); //Select automatic STMIN etc... calc. (AUTO 1) or manual (AUTO 0)
+ geant3->SetABAN(0); //Restore 3.16 behaviour for abandoned tracks
+ geant3->SetOPTI(2); //Select optimisation level for GEANT geometry searches (0,1,2)
+ geant3->SetERAN(5.e-7);
+
+ Float_t cut = 1.e-3; // 1MeV cut by default
+ Float_t tofmax = 1.e10;
+
+ // GAM ELEC NHAD CHAD MUON EBREM MUHAB EDEL MUDEL MUPA TOFMAX
+ geant3->SetCUTS(cut, cut, cut, cut, cut, cut, cut, cut, cut, cut,
+ tofmax);
+ //
+ //=======================================================================
+ // ************* STEERING parameters FOR ALICE SIMULATION **************
+ // --- Specify event type to be tracked through the ALICE setup
+ // --- All positions are in cm, angles in degrees, and P and E in GeV
+ if (gSystem->Getenv("CONFIG_NPARTICLES")){
+ int nParticles = atoi(gSystem->Getenv("CONFIG_NPARTICLES"));
+ }else{
+ int nParticles = 2;
+ } // end if
+//*********************************************
+// Example for Fixed Particle Gun
+//*********************************************
+ AliGenFixed *gener = new AliGenFixed(nParticles);
+ gener->SetMomentum(50);
+ gener->SetPhi(180.);
+ gener->SetTheta(95.);
+ gener->SetOrigin(0,0,0); //vertex position
+ gener->SetPart(13); //GEANT particle type
+ gener->Init();
+ //
+ // Activate this line if you want the vertex smearing to happen
+ // track by track
+ //
+ //gener->SetVertexSmear(perTrack);
+
+ gAlice->SetField(-999, 2); //Specify maximum magnetic field in Tesla (neg. ==> default field)
+
+ Int_t iABSO = 1;
+ Int_t iDIPO = 1;
+ Int_t iFMD = 1;
+ Int_t iFRAME = 1;
+ Int_t iHALL = 1;
+ Int_t iITS = 1;
+ Int_t iMAG = 1;
+ Int_t iMUON = 1;
+ Int_t iPHOS = 1;
+ Int_t iPIPE = 1;
+ Int_t iPMD = 1;
+ Int_t iRICH = 1;
+ Int_t iSHIL = 1;
+ Int_t iSTART = 1;
+ Int_t iTOF = 1;
+ Int_t iTPC = 1;
+ Int_t iTRD = 1;
+ Int_t iZDC = 1;
+ Int_t iEMCAL = 1;
+
+ //=================== Alice BODY parameters =============================
+ AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
+
+ if (iMAG){
+ //=================== MAG parameters ============================
+ // --- Start with Magnet since detector layouts may be depending ---
+ // --- on the selected Magnet dimensions ---
+ AliMAG *MAG = new AliMAG("MAG", "Magnet");
+ }
+ if (iABSO){
+ //=================== ABSO parameters ============================
+ AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber");
+ }
+ if (iDIPO){
+ //=================== DIPO parameters ============================
+ AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
+ }
+ if (iHALL){
+ //=================== HALL parameters ============================
+ AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
+ }
+ if (iFRAME){
+ //=================== FRAME parameters ============================
+ AliFRAME *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
+ }
+ if (iSHIL){
+ //=================== SHIL parameters ============================
+ AliSHIL *SHIL = new AliSHILvF("SHIL", "Shielding");
+ }
+ if (iPIPE){
+ //=================== PIPE parameters ============================
+ AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
+ }
+ if(iITS) {
+ //=================== ITS parameters ============================
+ //
+ // As the innermost detector in ALICE, the Inner Tracking System
+ // "impacts" on almost all other detectors. This involves the fact
+ // that the ITS geometry still has several options to be followed
+ //in parallel in order to determine the best set-up which minimizes
+ // the induced background. All the geometries available to date are
+ // described in the following. Read carefully the comments and use
+ // the default version (the only one uncommented) unless you are making
+ // comparisons and you know what you are doing. In this case just
+ // uncomment the ITS geometry you want to use and run Aliroot.
+ // Detailed geometries:
+ AliITSvPPRasymm *ITS = new AliITSvPPRasymm("ITS","New ITS PPR detailed version with asymmetric services");
+ ITS->SetMinorVersion(2); // don't touch this parameter if you're not an ITS developer
+ ITS->SetReadDet(kFALSE); // don't touch this parameter if you're not an ITS developer
+ ITS->SetThicknessDet1(200.); // detector thickness on layer 1 must be in the range [100,300]
+ ITS->SetThicknessDet2(200.); // detector thickness on layer 2 must be in the range [100,300]
+ ITS->SetThicknessChip1(200.); // chip thickness on layer 1 must be in the range [150,300]
+ ITS->SetThicknessChip2(200.);// chip thickness on layer 2 must be in the range [150,300]
+ ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ ITS->SetCoolingFluid(1); // 1 --> water ; 0 --> freon
+ // Geant3 <-> EUCLID conversion
+ // ============================
+ // SetEUCLID is a flag to output (=1) or not to output (=0) both
+ // geometry and media to two ASCII files (called by default
+ // ITSgeometry.euc and ITSgeometry.tme) in a format understandable
+ // to the CAD system EUCLID. The default (=0) means that you dont
+ // want to use this facility.
+ ITS->SetEUCLID(0);
+ }
+ if (iTPC){
+ //============================ TPC parameters =========================
+ // This allows the user to specify sectors for the SLOW (TPC geometry
+ // 2) Simulator. SecAL (SecAU) <0 means that ALL lower (upper)
+ // sectors are specified, any value other than that requires at least
+ // one sector (lower or upper)to be specified!
+ // Reminder: sectors 1-24 are lower sectors ( 1-12->z>0,13-24->z<0)
+ // sectors 25-72 are the upper ones (25-48->z>0,49-72->z<0)
+ // SecLows - number of lower sectors specified (up to 6)
+ // SecUps - number of upper sectors specified (up to 12)
+ // Sens - sensitive strips for the Slow Simulator !!!
+ // This does NOT work if all S or L-sectors are specified, i.e.
+ // if SecAL or SecAU < 0
+ //---------------------------------------------------------------------
+ // gROOT->LoadMacro("SetTPCParam.C");
+ // AliTPCParam *param = SetTPCParam();
+ AliTPC *TPC = new AliTPCv2("TPC", "Default");
+ // All sectors included
+ TPC->SetSecAL(-1);
+ TPC->SetSecAU(-1);
+ }
+ if (iTOF){
+ //=================== TOF parameters ============================
+ AliTOF *TOF = new AliTOFv2("TOF", "normal TOF");
+ }
+ if (iRICH){
+ //=================== RICH parameters ===========================
+ AliRICH *RICH = new AliRICHv3("RICH", "normal RICH");
+ }
+ if (iZDC){
+ //=================== ZDC parameters ============================
+ AliZDC *ZDC = new AliZDCv2("ZDC", "normal ZDC");
+ }
+ if (iTRD){
+ //=================== TRD parameters ============================
+ AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
+ // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 90% Xe +
+ // 10% CO2)
+ TRD->SetGasMix(1);
+ // With hole in front of PHOS
+ TRD->SetPHOShole();
+ // With hole in front of RICH
+ TRD->SetRICHhole();
+ // Switch on TR
+ AliTRDsim *TRDsim = TRD->CreateTR();
+ }
+ if (iFMD){
+ //=================== FMD parameters ============================
+ AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
+ FMD->SetRingsSi1(256);
+ FMD->SetRingsSi2(64);
+ FMD->SetSectorsSi1(20);
+ FMD->SetSectorsSi2(24);
+ }
+ if (iMUON){
+ //=================== MUON parameters ===========================
+ AliMUON *MUON = new AliMUONv1("MUON", "default");
+ }
+ if (iPHOS){
+ //=================== PHOS parameters ===========================
+ AliPHOS *PHOS = new AliPHOSv1("PHOS", "GPS2");
+ }
+ if (iPMD){
+ //=================== PMD parameters ============================
+ AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
+ PMD->SetPAR(1., 1., 0.8, 0.02);
+ PMD->SetIN(6., 18., -580., 27., 27.);
+ PMD->SetGEO(0.0, 0.2, 4.);
+ PMD->SetPadSize(0.8, 1.0, 1.0, 1.5);
+ }
+ if (iEMCAL){
+ //=================== EMCAL parameters ============================
+ AliEMCAL *EMCAL = new AliEMCALv1("EMCAL", "G56_2_55_19_104_14");
+ }
+ if (iSTART){
+ //=================== START parameters ============================
+ AliSTART *START = new AliSTARTv1("START", "START Detector");
+ }
+}
+//----------------------------------------------------------------------
+Float_t EtaToTheta(Float_t arg){
+ return (180./TMath::Pi())*2.*atan(exp(-arg));
+}
--- /dev/null
+void DrawS01()
+{
+ gMC->Gsatt("*", "seen", -1);
+ gMC->Gsatt("alic", "seen", 0);
+ gROOT->LoadMacro("ViewS01.C");
+ gInterpreter->ProcessLine("ViewS01()");
+ gMC->Gdopt("hide", "on");
+ gMC->Gdopt("shad", "on");
+ gMC->Gsatt("*", "fill", 7);
+ gMC->SetClipBox(".");
+ gMC->SetClipBox("*", 0, 100, -1000, 1000, -1000, 1000);
+ gMC->DefaultRange();
+ gMC->Gdraw("alic", 40, 30, 0, 10, 10, .1, .1);
+ gMC->Gdhead(1111, "Inner Tracking System");
+ gMC->Gdman(16, 6, "MAN");
+}
--- /dev/null
+void DrawSPD(TString view="Prospective" ){
+ AliITS *ITS = (AliITS*) gAlice->GetModule("ITS");
+ if(!ITS) {
+ cout << "You must initilize AliRoot with an ITS before plotting it"<<endl;
+ return;
+ } // end if
+ Int_t version1 = ITS->GetMajorVersion();
+ Int_t version2 = ITS->GetMinorVersion();
+ gMC->Gsatt("*", "seen", -1);
+ gMC->Gsatt("*", "fill", 7);
+ gMC->Gsatt("alic", "seen", 0);
+ gROOT->LoadMacro("ViewITSSPD.C");
+ ViewITSSPD(version1,version2);
+ //gInterpreter->ProcessLine("ViewSPD()");
+ gMC->Gdopt("hide", "on");
+ gMC->Gdopt("shad", "on");
+ gMC->SetClipBox(".");
+ if(view.Contains("Prospective")){
+ //SetClipBox(volume name,xmin,xmax,ymin,ymax,zmin,zmax)
+ gMC->SetClipBox("*",0.0,100.0,-1000.0,1000.0,0.0,1000.0);
+ gMC->DefaultRange();
+ //Gdraw(volume name,theta,phi,psi,axis originX,axis originY,scaleX,scaleY) angles in degrees
+ gMC->Gdraw("alic", 40.0, 30.0, 0.0, 11, 10, .4, .4);
+ gMC->Gdhead(1111, "Inner Tracking System");
+ //gMC->Gdman(16,6,"MAN");
+ } // end if
+ if(view.Contains("EndView")){
+ //SetClipBox(name,xmin,xmax,ymin,ymax,zmin,zmax)
+ gMC->SetClipBox("*",-100.0,100.0,-1000.0,1000.0,1.0,1000.0);
+ gMC->DefaultRange();
+ //Gdraw(volume name,theta,phi,psi,axis originX,axis originY,scaleX,scaleY) angles in degrees
+ gMC->Gdraw("alic",0.0,0.0,0.0,10.0,10.0,1.0,1.0);
+ gMC->Gdhead(1111, "Inner Tracking System");
+ //gMC->Gdman(16, 6, "MAN");
+ } // end if
+ if(view.Contains("SideView")){
+ //SetClipBox(name,xmin,xmax,ymin,ymax,zmin,zmax)
+ gMC->SetClipBox("*",0.0,100.0,-1000.0,1000.0,0.0,1000.0);
+ gMC->DefaultRange();
+ //Gdraw(volume name,theta,phi,psi,axis originX,axis originY,scaleX,scaleY) angles in degrees
+ gMC->Gdraw("alic",90.0,0.0,0.0,10.0,10.0,0.35,0.35);
+ gMC->Gdhead(1111, "Inner Tracking System");
+ //gMC->Gdman(16, 6, "MAN");
+ } // end if
+}
--- /dev/null
+////////////////////////////////////////////////////////////////////////
+//
+// name: MergeV1
+// date: 11.4.2002
+// last update: 11.4.2002
+// author: Jiri Chudoba
+// version: 1.0
+//
+// description:
+// creates digits from sdigits for several detectors
+// stores sdigits in separate file (or in the source file
+// with sdigits). Stores gAlice object and copies TE to the
+// file with digits
+// ITS region of Interest is set
+// test
+//
+// input:
+// TString fileNameSDigits ... input file with sdigits
+// TString fileNameDigits ... output file with digits
+// Int_t nEvents ... how many events to process
+// Int_t ITS, TPC, ... many flags for diff. detectors
+//
+// History:
+//
+// 04.04.02 - first version
+//
+////////////////////////////////////////////////////////////////////////
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "iostream.h"
+#include "TDatetime.h"
+#include "STEER/AliRun.h"
+#include "STEER/AliRunDigitizer.h"
+#include "ITS/AliITSDigitizer.h"
+#include "ITS/AliITS.h"
+#include "ITS/AliITSDetType.h"
+#include "ITS/AliITSresponseSDD.h"
+#include "TPC/AliTPCDigitizer.h"
+#include "TRD/AliTRDdigitizer.h"
+#include "PHOS/AliPHOSDigitizer.h"
+#include "MUON/AliMUONDigitizer.h"
+#include "RICH/AliRICHDigitizer.h"
+#include "TStopwatch.h"
+#endif
+
+// #include "AliHits2SDigits.C"
+
+// void AliCopyN(TString inputFile, TString outputFile);
+
+Int_t MergeV1(TString fileNameDigits="digits.root",
+ TString fileNameSDigitsSig="sig.sdigits.root",
+ TString fileNameSDigitsBgr="bgr.sdigits.root",
+ Int_t nEvents = 1, Int_t iITS = 2, Int_t iTPC = 0,
+ Int_t iTRD = 0, Int_t iPHOS = 0, Int_t iMUON = 0,
+ Int_t iRICH = 0, Int_t iCopy = 1)
+{
+// delete the current gAlice object, the one from input file
+// will be used
+
+ if(gAlice){
+ delete gAlice;
+ gAlice = 0;
+ } // end if gAlice
+
+ // Connect the Root Galice file containing Geometry, Kine and Hits
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(fileNameSDigitsSig.Data());
+ if(!file) file = new TFile(fileNameSDigitsSig.Data());
+ TDatime *ct0 = new TDatime(2002,04,26,00,00,00), ct = file->GetCreationDate();
+
+
+ // 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");
+ } // end if !gAlice
+
+ AliRunDigitizer * manager = new AliRunDigitizer(2,1);
+ manager->SetInputStream(0,fileNameSDigitsSig.Data());
+ manager->SetInputStream(1,fileNameSDigitsBgr.Data());
+ if (fileNameDigits != "") {
+// if (iCopy) {
+// AliCopyN(fileNameSDigitsSig,fileNameDigits);
+// }
+ manager->SetOutputFile(fileNameDigits);
+ }
+ manager->SetNrOfEventsToWrite(nEvents);
+
+ if (iITS) {
+ AliITSDigitizer *dITS = new AliITSDigitizer(manager);
+ if (iITS == 2) dITS->SetByRegionOfInterestFlag(1);
+ if(ct0->GetDate()>ct.GetDate()){
+ // For old files, must change SDD noise.
+ AliITS *ITS = (AliITS*) gAlice->GetDetector("ITS");
+ AliITSresponseSDD *resp1 = ITS->DetType(1)->GetResponseModel();
+ resp1->SetNoiseParam();
+ resp1->SetNoiseAfterElectronics();
+ Float_t n,b;
+ Int_t cPar[8];
+ resp1->GetNoiseParam(n,b);
+ n = resp1->GetNoiseAfterElectronics();
+ cPar[0]=0;
+ cPar[1]=0;
+ cPar[2]=(Int_t)(b + 2.*n + 0.5);
+ cPar[3]=(Int_t)(b + 2.*n + 0.5);
+ cPar[4]=0;
+ cPar[5]=0;
+ cPar[6]=0;
+ cPar[7]=0;
+ resp1->SetCompressParam(cPar);
+ } // end if
+ }
+ if (iTPC) AliTPCDigitizer *dTPC = new AliTPCDigitizer(manager);
+ if (iTRD) AliTRDdigitizer *dTRD = new AliTRDdigitizer(manager);
+ if (iPHOS) AliPHOSDigitizer *dPHOS = new AliPHOSDigitizer(manager);
+ if (iMUON) AliMUONDigitizer *dMUON = new AliMUONDigitizer(manager);
+ if (iRICH) AliRICHDigitizer *dRICH = new AliRICHDigitizer(manager);
+ TStopwatch timer;
+ timer.Start();
+ manager->Exec("deb all");
+ timer.Stop();
+ timer.Print();
+// delete gAlice;
+// gAlice = 0;
+ delete manager;
+}
+
+
+/*
+////////////////////////////////////////////////////////////////////////
+void AliCopyN(TString inputFileName, TString outputFileName) {
+// copy some objects
+
+ TFile *inputFile = OpenFile(inputFileName);
+ if (!inputFile) return;
+
+ TFile *outputFile = TFile::Open(outputFileName.Data(),"update");
+ if (!outputFile->IsOpen()) {
+ cerr<<"Can't open "<<outputFileName.Data()<<" !\n";
+ return;
+ }
+ if (!ImportgAlice(inputFile)) return;
+ AliCopy(inputFile, outputFile);
+ delete gAlice;
+ gAlice=0;
+ inputFile->Close();
+ delete inputFile;
+ outputFile->Close();
+ delete outputFile;
+}
+////////////////////////////////////////////////////////////////////////
+*/
--- /dev/null
+void ViewITSSPD(Int_t version1,Int_t version2){
+ gMC->Gsatt("ITSV","seen",0);// Air
+ gMC->Gsatt("ITSD","seen",0);// Air
+ gMC->Gsatt("IT12","seen",0);// Air
+ if(version1==10){gMC->Gsatt("I651","seen",1);gMC->Gsatt("I651","colo",7);} // Services
+ gMC->Gsatt("I12B","seen",0);// Air
+ gMC->Gsatt("I10B","seen",0);// Air
+ gMC->Gsatt("I105","seen",1);gMC->Gsatt("I105","colo",7);// SPD End ladder
+ gMC->Gsatt("I108","seen",1);gMC->Gsatt("I108","colo",7);// SPD bus
+ gMC->Gsatt("I109","seen",1);gMC->Gsatt("I109","colo",7);// SPD bus
+ gMC->Gsatt("I107","seen",0);//Air
+ gMC->Gsatt("I106","seen",1);gMC->Gsatt("I106","colo",6);// Silicon
+ gMC->Gsatt("I101","seen",1);gMC->Gsatt("I101","colo",6);// Silicon
+ gMC->Gsatt("ITS1","seen",1);gMC->Gsatt("ITS1","colo",6);// Silicon
+ gMC->Gsatt("I20B","seen",0);// Air
+ gMC->Gsatt("I105","seen",1);gMC->Gsatt("I105","colo",7);// SPD End ladder
+ gMC->Gsatt("I109","seen",1);gMC->Gsatt("I109","colo",7);// SPD bus
+ gMC->Gsatt("I108","seen",1);gMC->Gsatt("I108","colo",7);// SPD bus
+ gMC->Gsatt("I107","seen",0);// Air
+ gMC->Gsatt("I1D6","seen",1);gMC->Gsatt("I1D6","colo",6);// Silicon
+ gMC->Gsatt("I1D1","seen",1);gMC->Gsatt("I1D1","colo",6);// Silicon
+ gMC->Gsatt("ITS2","seen",1);gMC->Gsatt("ITS2","colo",6);// Silicon
+ gMC->Gsatt("I123","seen",1);gMC->Gsatt("I123","colo",1);// carbon
+ gMC->Gsatt("I121","seen",1);gMC->Gsatt("I121","colo",1);// carbon
+ gMC->Gsatt("I122","seen",1);gMC->Gsatt("I122","colo",1);// carbon
+ gMC->Gsatt("I120","seen",1);gMC->Gsatt("I120","colo",1);// carbon
+ gMC->Gsatt("I144","seen",1);gMC->Gsatt("I144","colo",1);// carbon
+ gMC->Gsatt("I143","seen",1);gMC->Gsatt("I143","colo",1);// carbon
+ gMC->Gsatt("I142","seen",1);gMC->Gsatt("I142","colo",1);// carbon
+ gMC->Gsatt("I141","seen",1);gMC->Gsatt("I141","colo",1);// carbon
+ gMC->Gsatt("I140","seen",1);gMC->Gsatt("I140","colo",1);// carbon
+ gMC->Gsatt("I139","seen",1);gMC->Gsatt("I139","colo",1);// carbon
+ gMC->Gsatt("I138","seen",1);gMC->Gsatt("I138","colo",1);// carbon
+ gMC->Gsatt("I137","seen",1);gMC->Gsatt("I137","colo",1);// carbon
+ gMC->Gsatt("I136","seen",1);gMC->Gsatt("I136","colo",1);// carbon
+ gMC->Gsatt("I135","seen",1);gMC->Gsatt("I135","colo",1);// carbon
+ gMC->Gsatt("I134","seen",1);gMC->Gsatt("I134","colo",1);// carbon
+ gMC->Gsatt("I133","seen",1);gMC->Gsatt("I133","colo",1);// carbon
+ gMC->Gsatt("I132","seen",1);gMC->Gsatt("I132","colo",1);// carbon
+ gMC->Gsatt("I131","seen",1);gMC->Gsatt("I131","colo",1);// carbon
+ gMC->Gsatt("I130","seen",1);gMC->Gsatt("I130","colo",1);// carbon
+ gMC->Gsatt("I129","seen",1);gMC->Gsatt("I129","colo",1);// carbon
+ gMC->Gsatt("I128","seen",1);gMC->Gsatt("I128","colo",1);// carbon
+ gMC->Gsatt("I126","seen",1);gMC->Gsatt("I126","colo",1);// carbon
+ gMC->Gsatt("I125","seen",1);gMC->Gsatt("I125","colo",1);// carbon
+ gMC->Gsatt("I124","seen",1);gMC->Gsatt("I124","colo",1);// carbon
+ gMC->Gsatt("I113","seen",0);//Air
+ gMC->Gsatt("I112","seen",1);gMC->Gsatt("I112","colo",1);// carbon
+ gMC->Gsatt("I111","seen",1);gMC->Gsatt("I111","colo",1);// carbon
+ gMC->Gsatt("I118","seen",1);gMC->Gsatt("I118","colo",2);// glue
+ gMC->Gsatt("I110","seen",1);gMC->Gsatt("I110","colo",1);// carbon
+ gMC->Gsatt("I114","seen",1);gMC->Gsatt("I114","colo",3);// INOX
+ gMC->Gsatt("I115","seen",1);gMC->Gsatt("I115","colo",4);// water
+ gMC->Gsatt("I116","seen",1);gMC->Gsatt("I116","colo",3);// INOX
+ gMC->Gsatt("I117","seen",1);gMC->Gsatt("I117","colo",4);// Water
+ gMC->Gsatt("I650","seen",0); //air
+ gMC->Gsatt("I666","seen",1);gMC->Gsatt("I666","colo",5);//Aluminum
+ gMC->Gsatt("I667","seen",1);gMC->Gsatt("I667","colo",5);//Aluminum
+ gMC->Gsatt("I668","seen",1);gMC->Gsatt("I668","colo",4);// water
+ gMC->Gsatt("I669","seen",1);gMC->Gsatt("I669","colo",3);// INOX
+ gMC->Gsatt("I670","seen",1);gMC->Gsatt("I670","colo",4);// water
+ gMC->Gsatt("I671","seen",1);gMC->Gsatt("I671","colo",6);// Silicon
+ gMC->Gsatt("I672","seen",1);gMC->Gsatt("I672","colo",4);// water
+ gMC->Gsatt("I673","seen",1);gMC->Gsatt("I673","colo",6);// Silicon
+ gMC->Gsatt("I674","seen",1);gMC->Gsatt("I674","colo",4);// water
+ gMC->Gsatt("I675","seen",1);gMC->Gsatt("I675","colo",5);//Aluminum
+ gMC->Gsatt("I676","seen",1);gMC->Gsatt("I676","colo",6);// Silicon
+ gMC->Gsatt("I677","seen",1);gMC->Gsatt("I677","colo",4);// water
+ gMC->Gsatt("I678","seen",1);gMC->Gsatt("I678","colo",5);//Aluminum
+}
--- /dev/null
+void ViewS01(){
+ gMC->Gsatt("ITSV","seen",0); // Air "Gen Air"
+ gMC->Gsatt("ITSD","seen",0); // Air "Gen Air"
+ gMC->Gsatt("IS01","seen",0); // Air "Air"
+ gMC->Gsatt("ICY2","seen",1); // Carbon Fiber "SDD C (M55J)
+ gMC->Gsatt("I212","seen",1); // Carbon "ITS SandW C"
+ gMC->Gsatt("I211","seen",1); // PC board "G10FR4"
+ gMC->Gsatt("I217","seen",1); // Carbon Fiber "SDD/SSD rings"
+ gMC->Gsatt("I219","seen",1); // Carbon Fiber "SDD/SSD rings"
+ gMC->Gsatt("I214","seen",1); // PC board "G10FR4"
+ gMC->Gsatt("I213","seen",1); // PC board "G10FR4"
+ gMC->Gsatt("I215","seen",2); // Air "Air"
+ gMC->Gsatt("I216","seen",1); // Carbon "ITS Sandw C"
+}
--- /dev/null
+void ViewSPD(){
+ gMC->Gsatt("ITSV","seen",0);
+ gMC->Gsatt("ITSD","seen",0);
+ gMC->Gsatt("IT12","seen",0);
+ gMC->Gsatt("I12B","seen",0);
+ gMC->Gsatt("I10B","seen",0);
+ gMC->Gsatt("I105","seen",1);
+ gMC->Gsatt("I108","seen",1);
+ gMC->Gsatt("I109","seen",1);
+ gMC->Gsatt("I107","seen",0);
+ gMC->Gsatt("I106","seen",1);
+ gMC->Gsatt("I101","seen",0);
+ gMC->Gsatt("ITS1","seen",1);
+ gMC->Gsatt("I20B","seen",0);
+ gMC->Gsatt("I105","seen",1);
+ gMC->Gsatt("I109","seen",1);
+ gMC->Gsatt("I108","seen",1);
+ gMC->Gsatt("I107","seen",0);
+ gMC->Gsatt("I1D6","seen",1);
+ gMC->Gsatt("I1D1","seen",0);
+ gMC->Gsatt("ITS2","seen",1);
+ gMC->Gsatt("I123","seen",1);
+ gMC->Gsatt("I121","seen",1);
+ gMC->Gsatt("I122","seen",1);
+ gMC->Gsatt("I120","seen",1);
+ gMC->Gsatt("I144","seen",1);
+ gMC->Gsatt("I143","seen",1);
+ gMC->Gsatt("I142","seen",1);
+ gMC->Gsatt("I141","seen",1);
+ gMC->Gsatt("I140","seen",1);
+ gMC->Gsatt("I139","seen",1);
+ gMC->Gsatt("I138","seen",1);
+ gMC->Gsatt("I137","seen",1);
+ gMC->Gsatt("I136","seen",1);
+ gMC->Gsatt("I135","seen",1);
+ gMC->Gsatt("I134","seen",1);
+ gMC->Gsatt("I133","seen",1);
+ gMC->Gsatt("I132","seen",1);
+ gMC->Gsatt("I131","seen",1);
+ gMC->Gsatt("I130","seen",1);
+ gMC->Gsatt("I129","seen",1);
+ gMC->Gsatt("I128","seen",1);
+ gMC->Gsatt("I126","seen",1);
+ gMC->Gsatt("I125","seen",1);
+ gMC->Gsatt("I124","seen",1);
+ gMC->Gsatt("I113","seen",0);
+ gMC->Gsatt("I112","seen",1);
+ gMC->Gsatt("I111","seen",1);
+ gMC->Gsatt("I118","seen",1);
+ gMC->Gsatt("I110","seen",1);
+ gMC->Gsatt("I114","seen",1);
+ gMC->Gsatt("I115","seen",1);
+ gMC->Gsatt("I116","seen",0);
+ gMC->Gsatt("I117","seen",1);
+}