[u/mrichter/AliRoot.git] / ITS / ITSgeoplot.C

#define __NOCOMPILED__
#ifdef __COMPILED__
#include<iostream.h>
#include<TROOT.h>
#include<TArrayI.h>
#include<TCanvas.h>
#include<TClassTable.h>
#include<TClonesArray.h>
#include<TFile.h>
#include<TH1.h>
#include<TH2.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 ITSgeoplot (char *opt="All+Rec", char *filename="galice.root") {
  /*******************************************************************
   *  This macro displays geometrical information related to the
   *  hits, digits and rec points in ITS.
   *  There are histograms that are not displayed (i.e. energy
   *  deposition) but that are saved onto a file (see below)
   *
   *  Options: Any combination of:
   *    1) subdetector name:  "SPD", "SDD", "SSD", "All" (default)
   *    2) Printouts:        "Verbose" Almost everything is printed out: 
   *                          it is wise to redirect the output onto a file
   *                    e.g.: .x ITSgeoplot.C("All+Verbose") > out.log 
   *    3) Rec Points option: "Rec"   ---> Uses Rec. Points (default)
   *                           otherwise ---> uses hits and digits only
   *    Examples:
   *       .x ITSgeoplot();  (All subdetectors; no-verbose; no-recpoints)
   *       .x ITSgeoplot("SPD+SSD+Verbose+Rec"); 
   *
   *  OUTPUT: It produces a root file with a list of histograms
   *          The list can be accessed either with the Root TBrowser
   *          or with the simple macro ITSgeoplotread.C
   *  WARNING: spatial information for SSD/DIGITS is obtained by pairing
   *          digits on p and n side originating from the same track, when
   *          possible. This (mis)use of DIGITS is tolerated for debugging 
   *          purposes only !!!!  The pairing in real life should be done
   *          starting from info really available... 
   * 
   *  COMPILATION: this macro can be compiled. 
   *      1)       Change the first line above from its default 
   *               value (#define __NOCOMPILED__) to
   *               #define __COMPILED__
   *      2)       You need to set your include path with
   * gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g");
   *      3)       If you are using root instead of aliroot you need to
   *               execute the macro loadlibs.C in advance
   *      4)       To compile this macro from root (or aliroot):
   *                 ---  .L ITSgeoplot.C++
   *                 ---  ITSgeoplot("ListOfParametersIfAny");
   *     
   *  M.Masera  14/05/2001 18:30
   *            Use DetToLocal instead of GetCxz for SPD and SDD
   ********************************************************************/

  extern void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb);
  extern Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
  extern void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos);
  //Options
  TString choice(opt);
  Bool_t All = choice.Contains("All");
  Bool_t verbose=choice.Contains("Verbose");
  Bool_t userec=choice.Contains("Rec");
  Int_t retcode=1; //return code
#ifdef __NOCOMPILED__
  if (gClassTable->GetID("AliRun") < 0) {
    gROOT->LoadMacro("loadlibs.C");
    loadlibs();
  }
  else {
#endif
    if(gAlice){
      delete gAlice;
      gAlice=0;
    }
#ifdef __NOCOMPILED__
  }
#endif
  // Connect the Root input  file containing Geometry, Kine and Hits
  // galice.root file by default

  TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
  if (!file) file = new TFile(filename);
  file->ls();

  // Get AliRun object from file

  if (!gAlice) {
    gAlice = (AliRun*)file->Get("gAlice");
    if (gAlice && verbose)cout<<"AliRun object found on file "<<filename<<endl;
    if(!gAlice){
      cout<<"Can't access AliRun object on file "<<filename<<endl;
      cout<<"Macro execution stopped!!!"<<endl;
      retcode=-1;
      return retcode;
    }
  }
  Int_t nparticles = gAlice->GetEvent(0);
  if(verbose) {
    cout<<" "<<endl<<" "<<endl;
    cout<<"******* Event processing started   *******"<<endl;
    cout<<"In the following, hits with 'StatusEntering' flag active"<<endl;
    cout<<"will not be processed"<<endl; 
    cout << "Number of particles=  " << nparticles <<endl;
  }

  // HITS
  TTree *TH = gAlice->TreeH();
  Stat_t ntracks = TH->GetEntries();
  if(verbose)cout<<"Number of primary tracks= "<<ntracks<<endl;

  // ITS
  AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
  Int_t nmodules;
  ITS->InitModules(-1,nmodules);
  cout<<"Number of ITS modules= "<<nmodules<<endl;
  cout<<"Filling modules... It takes a while, now. Please be patient"<<endl;
  ITS->FillModules(0,0,nmodules," "," ");
  cout<<"ITS modules .... DONE!"<<endl;

  // DIGITS
  TTree *TD = gAlice->TreeD();

  //RECPOINTS
  TTree *TR = gAlice->TreeR();
  TClonesArray *ITSrec  = ITS->RecPoints();
  if(userec && (!TR || !ITSrec)){
    userec = kFALSE;
    cout<<" "<<endl<<"======================================================="<<endl;
    cout<<"WARNING: there are no RECPOINTS on this file ! "<<endl;
    cout<<"======================================================="<<endl<<" "<<endl;
  }
  //local variables
  Int_t mod;   //module number
  Int_t nbytes = 0; 
  Double_t pos[3];  // Global position of the current module
  Float_t ragdet; // Radius of detector (x y plane)
  Int_t first; // first module
  Int_t last; // last module
  Int_t nrecp; //number of RecPoints for a given module

  //List of histograms
  TObjArray histos(26,0);  // contains the pointers to the histograms
  // Book histograms SPD
  TH2F *rspd = new TH2F("rspd","Radii of digits - SPD",50,-10.,10.,50,-10.,10.);
  TH2F *rhspd = new TH2F("rhspd","Radii of hits - SPD",50,-10.,10.,50,-10.,10.);
  TH2F *rmspd = new TH2F("rmspd","Radii of SPD modules",50,-10.,10.,50,-10.,10.);
  TH1F *zspd = new TH1F("zspd","Z of digits - SPD",100,-30.,30.);
  TH1F *zhspd = new TH1F("zhspd","Z of hits - SPD",100,-30.,30.);
  TH1F *zmspd = new TH1F("zmspd","Z of SPD modules",100,-30,30.);
  TH2F *rrspd = new TH2F("rrspd","Radii of recpoints - SPD",50,-10.,10.,50,-10.,10.);
  TH1F *zrspd = new TH1F("zrspd","Z of recpoints - SPD",100,-30.,30.);
  TH1F *enespd = new TH1F("enespd","Energy deposition SPD (KeV)",100,0.,1000.);
  histos.AddLast(rspd);  // 0
  histos.AddLast(rhspd); // 1
  histos.AddLast(rmspd); // 2
  histos.AddLast(zspd);  // 3
  histos.AddLast(zhspd); // 4
  histos.AddLast(zmspd); // 5
  histos.AddLast(rrspd); // 6
  histos.AddLast(zrspd); // 7
  histos.AddLast(enespd); // 8
  // Book histograms SDD
  TH2F *rsdd = new TH2F("rsdd","Radii of digits - SDD",50,-40.,40.,50,-40.,40.);
  TH2F *rhsdd = new TH2F("rhsdd","Radii of hits - SDD",50,-40.,40.,50,-40.,40.);
  TH2F *rmsdd = new TH2F("rmsdd","Radii of SDD modules",50,-40.,40.,50,-40.,40.);
  TH1F *zsdd = new TH1F("zsdd","Z of digits - SDD",100,-40.,40.);
  TH1F *zhsdd = new TH1F("zhsdd","Z of hits - SDD",100,-40.,40.);
  TH1F *zmsdd = new TH1F("zmsdd","Z of SDD modules",100,-40,40.);
  TH2F *rrsdd = new TH2F("rrsdd","Radii of recpoints - SDD",50,-40.,40.,50,-40.,40.);
  TH1F *zrsdd = new TH1F("zrsdd","Z of recpoints - SDD",100,-40.,40.);
  TH1F *enesdd = new TH1F("enesdd","Energy deposition SDD (KeV)",100,0.,1000.);
  histos.AddLast(rsdd);  // 9
  histos.AddLast(rhsdd); // 10
  histos.AddLast(rmsdd); // 11
  histos.AddLast(zsdd);  // 12
  histos.AddLast(zhsdd); // 13
  histos.AddLast(zmsdd); // 14
  histos.AddLast(rrsdd); // 15
  histos.AddLast(zrsdd); // 16
  histos.AddLast(enesdd); // 17
  // Book histogram SSD
  TH2F *rssd = new TH2F("rssd","Radii of digits - SSD",50,-50.,50.,50,-50.,50.);
  TH2F *rhssd = new TH2F("rhssd","Radii of hits - SSD",50,-50.,50.,50,-50.,50.);
  TH2F *rmssd = new TH2F("rmssd","Radii of SSD modules",50,-50.,50.,50,-50.,50.);
  TH1F *zssd = new TH1F("zssd","Z of digits - SSD",100,-70.,70.);
  TH1F *zhssd = new TH1F("zhssd","Z of hits - SSD",100,-70.,70.);
  TH1F *zmssd = new TH1F("zmssd","Z of SSD modules",100,-70,70.);
  TH2F *rrssd = new TH2F("rrssd","Radii of recpoints - SSD",50,-50.,50.,50,-50.,50.);
  TH1F *zrssd = new TH1F("zrssd","Z of recpoints - SSD",100,-70.,70.);
  TH1F *enessd = new TH1F("enessd","Energy deposition SSD (KeV)",100,0.,1000.);
  histos.AddLast(rssd);  // 18
  histos.AddLast(rhssd); // 19
  histos.AddLast(rmssd); // 20
  histos.AddLast(zssd);  // 21
  histos.AddLast(zhssd); // 22
  histos.AddLast(zmssd); // 23
  histos.AddLast(rrssd); // 24
  histos.AddLast(zrssd); // 25
  histos.AddLast(enessd); // 26
  //
  // Loop on subdetectors
  // 
  AliITSgeom *geom = ITS->GetITSgeom();
  TString detna; // subdetector name
  for(Int_t subd=0;subd<3;subd++){
    if(All || (choice.Contains("SPD") && subd==0) || (choice.Contains("SDD") && subd==1) || (choice.Contains("SSD") && subd==2)){
      // Prepare array for the digits
      TClonesArray *ITSdigits  = ITS->DigitsAddress(subd);
      Bool_t usedigits = kTRUE;
      if(!ITSdigits){
        usedigits = kFALSE;
        cout<<" "<<endl<<"======================================================="<<endl;
        cout<<"WARNING: there are no DIGITS on this file ! "<<endl;
        cout<<"======================================================="<<endl<<" "<<endl;
      }
      // Get segmentation model
      AliITSDetType *iDetType=ITS->DetType(subd);
      if(subd==0)detna="SPD";
      if(subd==1)detna="SDD";
      if(subd==2)detna="SSD";
      AliITSsegmentation *seg=(AliITSsegmentation*)iDetType->GetSegmentationModel();
      // Loop on modules
      first = geom->GetStartDet(subd);
      last = geom->GetLastDet(subd);
      if(verbose){
        cout<<"     "<<endl<<"-------------------------------------"<<endl;
        cout<<"Start processing subdetector "<<detna<<endl;
        cout<<detna<<" first module "<<first<<endl;
        cout<<detna<<" last module "<<last<<endl;
        cout<<" "<<endl<<" "<<endl;
      }
      for (mod=first; mod<=last; mod++){
        geom->GetTrans(mod,pos);  // position of the module in the MRS
        ragdet=sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
        // The following 2 histos are a check of the geometry
        TH2F *bidi = (TH2F*)histos.At(2+subd*9);
        TH1F *uni = (TH1F*)histos.At(5+subd*9);
        bidi->Fill(pos[0],pos[1]);
        uni->Fill(pos[2]);
        if(verbose){
          cout<<"==========================================================="<<endl;
          cout<<detna<<" module="<<mod<<endl;
          cout<<"Mod. coordinates: "<<pos[0]<<", "<<pos[1]<<", "<<pos[2]<<" Radius "<<ragdet<<endl;
        }

        // Hits
        GetHitsCoor(ITS,mod,histos,subd,verbose);

        //RecPoints     
        if(userec){
          ITS->ResetRecPoints();
          TR->GetEvent(mod);
          TH2F *bidi=(TH2F*)histos.At(6+subd*9);
          TH1F *uni=(TH1F*)histos.At(7+subd*9);
          nrecp=GetRecCoor(geom,ITSrec,mod,bidi,uni,verbose);
        }
     
        // Digits
        if(usedigits){
          ITS->ResetDigits();
          nbytes += TD->GetEvent(mod);
          GetDigits(seg,geom,ITSdigits,subd,mod,verbose,histos);
        }

      } // End of loop on the modules
      TH1F *h1tmp;
      TH2F *h2tmp;
      //  Plot the histograms
      TCanvas *current=0; // current Canvas (1--> SPD, 2---> SDD, 3---> SSD)
      if(subd==0){
        // Prepare canvas 1
        TCanvas *c1 = new TCanvas("c1","SPD",10,10,600,900);
        c1->Divide(2,3);
        current=c1;
      }
      if(subd==1){
        // Prepare canvas 2
        TCanvas *c2 = new TCanvas("c2","SDD",40,40,600,900);
        c2->Divide(2,3);
        current=c2;
      }
      if(subd==2){
        // Prepare canvas 3
        TCanvas *c3 = new TCanvas("c3","SSD",70,70,600,900);
        c3->Divide(2,3);
        current=c3;
      }
      current->cd(1);
      h2tmp = (TH2F*)histos.At(9*subd);
      h2tmp->Draw();
      current->cd(2);
      h1tmp=(TH1F*)histos.At(3+subd*9);
      h1tmp->Draw();
      current->cd(3);
      h2tmp=(TH2F*)histos.At(1+9*subd);
      h2tmp->Draw();
      current->cd(4);
      h1tmp=(TH1F*)histos.At(4+subd*9);
      h1tmp->Draw();
   
      if(userec){
        current->cd(5);
        h2tmp=(TH2F*)histos.At(6+9*subd);
        h2tmp->Draw();
        current->cd(6);
        h1tmp=(TH1F*)histos.At(7+subd*9);
        h1tmp->Draw();
      }
  
      else {
        current->cd(5);
        h2tmp=(TH2F*)histos.At(2+9*subd);
        h2tmp->Draw();
        current->cd(6);
        h2tmp=(TH2F*)histos.At(5+9*subd);
        h2tmp->Draw();
      }
    } // if(All.....
  } // end of loop on subdetectors
  // Save the histograms
  TFile *fh = new TFile("ITSgeoplot.root","recreate");
  // The list is written to file as a single entry
  TList *lihi = new TList();
  // copy the pointers to the histograms to a TList object.
  // The histograms concerning recpoints are not copied if
  // 'userec' is false.
  for(Int_t i=0;i<histos.GetEntriesFast();i++){
    if(choice.Contains("All") || (choice.Contains("SPD") && i<8) || (choice.Contains("SDD") && i>7 && i<16) || (choice.Contains("SSD") && i>15)){
      if(!(!userec && ((i+2)%9==0 || (i+1)%9==0)))lihi->Add(histos.At(i));
    }
  }
  lihi->Write("Histograms ITS hits+digits+recpoints",TObject::kSingleKey);
  fh->Close();

  return retcode;
}


void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb){
  TH2F *h2=(TH2F*)histos.At(1+subd*9);
  TH1F *h1=(TH1F*)histos.At(4+subd*9);
  TH1F *ener=(TH1F*)histos.At(8+subd*9);
  AliITS *ITS= (AliITS*)its;
  AliITSmodule *modu = ITS->GetModule(mod);
  TObjArray *fHits = modu->GetHits();
  Int_t nhits = fHits->GetEntriesFast();
  if(nhits>0){
    if(verb){
      cout<<"-------------------------------------------------------"<<endl;
      cout<<"Number of HITS for module "<<mod<<": "<<nhits<<endl;
    }
    for (Int_t hit=0;hit<nhits;hit++) {
      AliITShit *iHit = (AliITShit*) fHits->At(hit);
      if(!iHit->StatusEntering()){
        Float_t x=iHit->GetXG();
        Float_t y=iHit->GetYG();
        Float_t z=iHit->GetZG();
        Float_t edep=iHit->GetIonization()*1000000;
        h2->Fill(x,y);
        h1->Fill(z);
        ener->Fill(edep);
        if(verb){
          cout<<"hit # "<<hit<<" Global coordinates "<<x<<" "<<y<<" "<<z<<endl;
          cout<<"track # "<<iHit->GetTrack()<<" energy deposition (KeV)= "<<edep<<endl;
        }
      }
    }
  }
}


Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
  AliITSgeom *geom = (AliITSgeom*)ge;
  Int_t nrecp = ITSrec->GetEntries();
  if(nrecp>0){
    Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
    Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
    if(verb){
      cout<<"-------------------------------------------------------"<<endl;
      cout<<"Number of REC POINTS for module "<<mod<<": "<<nrecp<<endl;
    }
    for(Int_t irec=0;irec<nrecp;irec++) {
      AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
      lc[0]=recp->GetX();
      lc[2]=recp->GetZ();
      geom->LtoG(mod,lc,gc);
      if(verb){
        cout<<"recp # "<<irec<<" local coordinates. lx= "<<lc[0]<<" lz= "<<lc[2]<<endl;
        Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
        cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<r<<endl;
      }
      h2->Fill(gc[0],gc[1]);
      h1->Fill(gc[2]);
    }
  }
  return nrecp;
}

void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos){
  AliITSsegmentation *seg = (AliITSsegmentation*)tmps;
  AliITSgeom *geom = (AliITSgeom*)ge;
  Int_t layer;
  Int_t ladder;
  Int_t detec;
  if(subd==2){
    geom->GetModuleId(mod,layer,ladder,detec);
    seg->SetLayer(layer);
  }
  Float_t lcoor[3]; for(Int_t j=0; j<3; j++) lcoor[j]=0.;  //local coord dig.
  Float_t gcoor[3]; for(Int_t j=0; j<3; j++) gcoor[j]=0.; // global coo. dig.
  Float_t ragdig; // Radius digit
  TArrayI ssdone(5000);  // used to match p and n side digits of strips
  TArrayI pair(5000);    // as above 
  Int_t ndigits = ITSdigits->GetEntries();
  AliITSdigit *digs;
  if(ndigits){ 
    if(verbose){
      cout<<"-------------------------------------------------------"<<endl;
      cout<<"Number of DIGITS for module "<<mod<<": "<<ndigits<<endl;
    }
    // Get the coordinates of the module
    if(subd==2){
      for (Int_t digit=0;digit<ndigits;digit++){
        ssdone[digit]=0;
        pair[digit]=0;
      }
    }
    for (Int_t digit=0;digit<ndigits;digit++) {
      digs = (AliITSdigit*)ITSdigits->UncheckedAt(digit);
      Int_t iz=digs->fCoord1;  // cell number z
      Int_t ix=digs->fCoord2;  // cell number x
      // Get local coordinates of the element 
      if(subd<2){
        seg->DetToLocal(ix,iz,lcoor[0],lcoor[2]);
      }
      else{
        // SSD: if iz==0 ---> N side; if iz==1 P side
        if(ssdone[digit]==0){
          ssdone[digit]=1;
          pair[digit]=-1;
          Bool_t pside=(iz==1);
          Bool_t impaired=kTRUE;
          Int_t pstrip=0;
          Int_t nstrip=0;
          if(pside)pstrip=ix;
          if(!pside)nstrip=ix;
          for(Int_t digi2=0;digi2<ndigits;digi2++){
            if(ssdone[digi2]==0 && impaired){
              AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(digi2);
              if(dig2->fCoord1 != iz && dig2->GetTrack(0)==digs->GetTrack(0) && dig2->GetTrack(0)>0){
                ssdone[digi2]=2;
                pair[digit]=digi2;
                if(pside)nstrip=dig2->fCoord2;
                if(!pside)pstrip=dig2->fCoord2;
                impaired=kFALSE;
              }
            }
          }
          if(!impaired)seg->GetPadCxz(pstrip,nstrip,lcoor[0],lcoor[2]);
        }
      }
      if(subd<2 || (subd==2 && ssdone[digit]==1)){
        Int_t coor1=digs->fCoord1;
        Int_t coor2=digs->fCoord2;
        Int_t tra0=digs->GetTrack(0);
        if(verbose){
          cout<<"digit # "<<digit<<" fCoord1= "<<coor1<<" fCoord2= "<<coor2<<" track "<<tra0<<" "<<endl;
          if(subd<2)cout<<"local coordinates -- x="<<lcoor[0]<<", z="<<lcoor[2]<<endl;
          if(subd==2){
            if(pair[digit]==-1){
              cout<<"(digit  not paired)"<<endl;
            }
            else {
              cout<<"local coordinates -- x="<<lcoor[0]<<", z="<<lcoor[2]<<endl;
              Int_t dtmp=pair[digit];
              AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(dtmp);
              Int_t coor1b=dig2->fCoord1;
              Int_t coor2b=dig2->fCoord2;
              Int_t tra0b=dig2->GetTrack(0);
              cout<<"(digit paired with digit #"<<dtmp<<endl;
              cout<<"with fCoord1= "<<coor1b<<" fCoord2= "<<coor2b<<" track "<<tra0b<<")"<<endl;
            }
          }
        }
        if(subd<2 || (subd==2 && pair[digit]!=-1)){
          // Global coordinates of the element
          //SDD and SPD use cm, SSD microns (GetPadCxz)
          if(subd==2)for(Int_t j=0;j<3;j++)lcoor[j]=lcoor[j]/10000.;
          lcoor[1]=0.;
          geom->LtoG(mod,lcoor,gcoor);  // global coord. in cm
          ragdig=sqrt(gcoor[0]*gcoor[0]+gcoor[1]*gcoor[1]);
          if(verbose)cout<<"global coordinates "<<gcoor[0]<<" "<<gcoor[1]<<" "<<gcoor[2]<<" Radius "<<ragdig<<endl;
          //Fill histograms
          TH2F *bidi = (TH2F*)histos.At(subd*9);
          TH1F *uni = (TH1F*)histos.At(3+subd*9);
          bidi->Fill(gcoor[0],gcoor[1]);
          uni->Fill(gcoor[2]);
        }
      }
    } // loop on digits for this module
  } // if(ndigits>0....
}
Commit	Line	Data
79f97fd2	1	#define __NOCOMPILED__
	2	#ifdef __COMPILED__
	3	#include<iostream.h>
	4	#include<TROOT.h>
	5	#include<TArrayI.h>
	6	#include<TCanvas.h>
	7	#include<TClassTable.h>
	8	#include<TClonesArray.h>
	9	#include<TFile.h>
	10	#include<TH1.h>
	11	#include<TH2.h>
	12	#include<TObject.h>
	13	#include<TObjArray.h>
	14	#include<TTree.h>
	15	#include <AliRun.h>
	16	#include <AliITS.h>
	17	#include <AliITSgeom.h>
	18	#include <AliITSDetType.h>
	19	#include <AliITSRecPoint.h>
	20	#include <AliITSdigit.h>
	21	#include <AliITShit.h>
	22	#include <AliITSmodule.h>
	23	#include <AliITSsegmentation.h>
	24	#include <AliITSsegmentationSPD.h>
	25	#include <AliITSsegmentationSDD.h>
	26	#include <AliITSsegmentationSSD.h>
	27	#endif
2f101c6a	28	Int_t ITSgeoplot (char opt="All+Rec", char filename="galice.root") {
	29	/*******************************************************************
	30	* This macro displays geometrical information related to the
	31	* hits, digits and rec points in ITS.
	32	* There are histograms that are not displayed (i.e. energy
	33	* deposition) but that are saved onto a file (see below)
	34	*
	35	* Options: Any combination of:
	36	* 1) subdetector name: "SPD", "SDD", "SSD", "All" (default)
	37	* 2) Printouts: "Verbose" Almost everything is printed out:
	38	* it is wise to redirect the output onto a file
	39	* e.g.: .x ITSgeoplot.C("All+Verbose") > out.log
	40	* 3) Rec Points option: "Rec" ---> Uses Rec. Points (default)
	41	* otherwise ---> uses hits and digits only
	42	* Examples:
	43	* .x ITSgeoplot(); (All subdetectors; no-verbose; no-recpoints)
	44	* .x ITSgeoplot("SPD+SSD+Verbose+Rec");
	45	*
	46	* OUTPUT: It produces a root file with a list of histograms
	47	* The list can be accessed either with the Root TBrowser
	48	* or with the simple macro ITSgeoplotread.C
	49	* WARNING: spatial information for SSD/DIGITS is obtained by pairing
	50	* digits on p and n side originating from the same track, when
	51	* possible. This (mis)use of DIGITS is tolerated for debugging
	52	* purposes only !!!! The pairing in real life should be done
79f97fd2	53	* starting from info really available...
	54	*
	55	* COMPILATION: this macro can be compiled.
	56	* 1) Change the first line above from its default
	57	* value (#define __NOCOMPILED__) to
	58	* #define __COMPILED__
	59	* 2) You need to set your include path with
	60	* gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g");
	61	* 3) If you are using root instead of aliroot you need to
	62	* execute the macro loadlibs.C in advance
	63	* 4) To compile this macro from root (or aliroot):
	64	* --- .L ITSgeoplot.C++
	65	* --- ITSgeoplot("ListOfParametersIfAny");
2f101c6a	66	*
faafc0eb	67	* M.Masera 14/05/2001 18:30
faafc0eb	68	* Use DetToLocal instead of GetCxz for SPD and SDD
2f101c6a	69	********************************************************************/
	70
	71	extern void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb);
	72	extern Int_t GetRecCoor(TObject ge, TClonesArray ITSrec, Int_t mod, TH2F h2, TH1F h1, Bool_t verb);
	73	extern void GetDigits(TObject tmps,TObject ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos);
2f101c6a	74	//Options
	75	TString choice(opt);
	76	Bool_t All = choice.Contains("All");
	77	Bool_t verbose=choice.Contains("Verbose");
	78	Bool_t userec=choice.Contains("Rec");
	79	Int_t retcode=1; //return code
79f97fd2	80	#ifdef __NOCOMPILED__
2f101c6a	81	if (gClassTable->GetID("AliRun") < 0) {
	82	gROOT->LoadMacro("loadlibs.C");
	83	loadlibs();
	84	}
	85	else {
79f97fd2	86	#endif
2f101c6a	87	if(gAlice){
	88	delete gAlice;
	89	gAlice=0;
	90	}
79f97fd2	91	#ifdef __NOCOMPILED__
2f101c6a	92	}
79f97fd2	93	#endif
2f101c6a	94	// Connect the Root input file containing Geometry, Kine and Hits
	95	// galice.root file by default
	96
	97	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject(filename);
	98	if (!file) file = new TFile(filename);
	99	file->ls();
	100
	101	// Get AliRun object from file
	102
	103	if (!gAlice) {
	104	gAlice = (AliRun*)file->Get("gAlice");
	105	if (gAlice && verbose)cout<<"AliRun object found on file "<<filename<<endl;
	106	if(!gAlice){
	107	cout<<"Can't access AliRun object on file "<<filename<<endl;
	108	cout<<"Macro execution stopped!!!"<<endl;
	109	retcode=-1;
	110	return retcode;
	111	}
	112	}
	113	Int_t nparticles = gAlice->GetEvent(0);
	114	if(verbose) {
	115	cout<<" "<<endl<<" "<<endl;
	116	cout<<"***** Event processing started *****"<<endl;
	117	cout<<"In the following, hits with 'StatusEntering' flag active"<<endl;
	118	cout<<"will not be processed"<<endl;
	119	cout << "Number of particles= " << nparticles <<endl;
	120	}
	121
	122	// HITS
	123	TTree *TH = gAlice->TreeH();
	124	Stat_t ntracks = TH->GetEntries();
	125	if(verbose)cout<<"Number of primary tracks= "<<ntracks<<endl;
	126
	127	// ITS
	128	AliITS ITS = (AliITS)gAlice->GetModule("ITS");
	129	Int_t nmodules;
	130	ITS->InitModules(-1,nmodules);
	131	cout<<"Number of ITS modules= "<<nmodules<<endl;
	132	cout<<"Filling modules... It takes a while, now. Please be patient"<<endl;
	133	ITS->FillModules(0,0,nmodules," "," ");
	134	cout<<"ITS modules .... DONE!"<<endl;
	135
	136	// DIGITS
	137	TTree *TD = gAlice->TreeD();
	138
	139	//RECPOINTS
	140	TTree *TR = gAlice->TreeR();
	141	TClonesArray *ITSrec = ITS->RecPoints();
	142	if(userec && (!TR \|\| !ITSrec)){
	143	userec = kFALSE;
	144	cout<<" "<<endl<<"======================================================="<<endl;
	145	cout<<"WARNING: there are no RECPOINTS on this file ! "<<endl;
	146	cout<<"======================================================="<<endl<<" "<<endl;
	147	}
	148	//local variables
	149	Int_t mod; //module number
	150	Int_t nbytes = 0;
	151	Double_t pos[3]; // Global position of the current module
	152	Float_t ragdet; // Radius of detector (x y plane)
	153	Int_t first; // first module
	154	Int_t last; // last module
	155	Int_t nrecp; //number of RecPoints for a given module
	156
	157	//List of histograms
158	TObjArray histos(26,0); // contains the pointers to the histograms
159	// Book histograms SPD
160	TH2F *rspd = new TH2F("rspd","Radii of digits - SPD",50,-10.,10.,50,-10.,10.);
161	TH2F *rhspd = new TH2F("rhspd","Radii of hits - SPD",50,-10.,10.,50,-10.,10.);
162	TH2F *rmspd = new TH2F("rmspd","Radii of SPD modules",50,-10.,10.,50,-10.,10.);
163	TH1F *zspd = new TH1F("zspd","Z of digits - SPD",100,-30.,30.);
164	TH1F *zhspd = new TH1F("zhspd","Z of hits - SPD",100,-30.,30.);
165	TH1F *zmspd = new TH1F("zmspd","Z of SPD modules",100,-30,30.);
166	TH2F *rrspd = new TH2F("rrspd","Radii of recpoints - SPD",50,-10.,10.,50,-10.,10.);
167	TH1F *zrspd = new TH1F("zrspd","Z of recpoints - SPD",100,-30.,30.);
168	TH1F *enespd = new TH1F("enespd","Energy deposition SPD (KeV)",100,0.,1000.);
169	histos.AddLast(rspd); // 0
170	histos.AddLast(rhspd); // 1
171	histos.AddLast(rmspd); // 2
172	histos.AddLast(zspd); // 3
173	histos.AddLast(zhspd); // 4
174	histos.AddLast(zmspd); // 5
175	histos.AddLast(rrspd); // 6
176	histos.AddLast(zrspd); // 7
177	histos.AddLast(enespd); // 8
178	// Book histograms SDD
179	TH2F *rsdd = new TH2F("rsdd","Radii of digits - SDD",50,-40.,40.,50,-40.,40.);
180	TH2F *rhsdd = new TH2F("rhsdd","Radii of hits - SDD",50,-40.,40.,50,-40.,40.);
181	TH2F *rmsdd = new TH2F("rmsdd","Radii of SDD modules",50,-40.,40.,50,-40.,40.);
182	TH1F *zsdd = new TH1F("zsdd","Z of digits - SDD",100,-40.,40.);
183	TH1F *zhsdd = new TH1F("zhsdd","Z of hits - SDD",100,-40.,40.);
184	TH1F *zmsdd = new TH1F("zmsdd","Z of SDD modules",100,-40,40.);
185	TH2F *rrsdd = new TH2F("rrsdd","Radii of recpoints - SDD",50,-40.,40.,50,-40.,40.);
186	TH1F *zrsdd = new TH1F("zrsdd","Z of recpoints - SDD",100,-40.,40.);
187	TH1F *enesdd = new TH1F("enesdd","Energy deposition SDD (KeV)",100,0.,1000.);
188	histos.AddLast(rsdd); // 9
189	histos.AddLast(rhsdd); // 10
190	histos.AddLast(rmsdd); // 11
191	histos.AddLast(zsdd); // 12
192	histos.AddLast(zhsdd); // 13
193	histos.AddLast(zmsdd); // 14
194	histos.AddLast(rrsdd); // 15
195	histos.AddLast(zrsdd); // 16
196	histos.AddLast(enesdd); // 17
197	// Book histogram SSD
198	TH2F *rssd = new TH2F("rssd","Radii of digits - SSD",50,-50.,50.,50,-50.,50.);
199	TH2F *rhssd = new TH2F("rhssd","Radii of hits - SSD",50,-50.,50.,50,-50.,50.);
200	TH2F *rmssd = new TH2F("rmssd","Radii of SSD modules",50,-50.,50.,50,-50.,50.);
201	TH1F *zssd = new TH1F("zssd","Z of digits - SSD",100,-70.,70.);
202	TH1F *zhssd = new TH1F("zhssd","Z of hits - SSD",100,-70.,70.);
203	TH1F *zmssd = new TH1F("zmssd","Z of SSD modules",100,-70,70.);
204	TH2F *rrssd = new TH2F("rrssd","Radii of recpoints - SSD",50,-50.,50.,50,-50.,50.);
205	TH1F *zrssd = new TH1F("zrssd","Z of recpoints - SSD",100,-70.,70.);
206	TH1F *enessd = new TH1F("enessd","Energy deposition SSD (KeV)",100,0.,1000.);
207	histos.AddLast(rssd); // 18
208	histos.AddLast(rhssd); // 19
209	histos.AddLast(rmssd); // 20
210	histos.AddLast(zssd); // 21
211	histos.AddLast(zhssd); // 22
212	histos.AddLast(zmssd); // 23
213	histos.AddLast(rrssd); // 24
214	histos.AddLast(zrssd); // 25
215	histos.AddLast(enessd); // 26
216	//
217	// Loop on subdetectors
218	//
219	AliITSgeom *geom = ITS->GetITSgeom();
79f97fd2	220	TString detna; // subdetector name
2f101c6a	221	for(Int_t subd=0;subd<3;subd++){
	222	if(All \|\| (choice.Contains("SPD") && subd==0) \|\| (choice.Contains("SDD") && subd==1) \|\| (choice.Contains("SSD") && subd==2)){
	223	// Prepare array for the digits
	224	TClonesArray *ITSdigits = ITS->DigitsAddress(subd);
	225	Bool_t usedigits = kTRUE;
	226	if(!ITSdigits){
	227	usedigits = kFALSE;
	228	cout<<" "<<endl<<"======================================================="<<endl;
	229	cout<<"WARNING: there are no DIGITS on this file ! "<<endl;
	230	cout<<"======================================================="<<endl<<" "<<endl;
	231	}
	232	// Get segmentation model
	233	AliITSDetType *iDetType=ITS->DetType(subd);
2f101c6a	234	if(subd==0)detna="SPD";
	235	if(subd==1)detna="SDD";
	236	if(subd==2)detna="SSD";
	237	AliITSsegmentation seg=(AliITSsegmentation)iDetType->GetSegmentationModel();
	238	// Loop on modules
	239	first = geom->GetStartDet(subd);
	240	last = geom->GetLastDet(subd);
	241	if(verbose){
	242	cout<<" "<<endl<<"-------------------------------------"<<endl;
	243	cout<<"Start processing subdetector "<<detna<<endl;
	244	cout<<detna<<" first module "<<first<<endl;
	245	cout<<detna<<" last module "<<last<<endl;
	246	cout<<" "<<endl<<" "<<endl;
	247	}
	248	for (mod=first; mod<=last; mod++){
	249	geom->GetTrans(mod,pos); // position of the module in the MRS
	250	ragdet=sqrt(pos[0]pos[0]+pos[1]pos[1]);
	251	// The following 2 histos are a check of the geometry
	252	TH2F bidi = (TH2F)histos.At(2+subd*9);
	253	TH1F uni = (TH1F)histos.At(5+subd*9);
	254	bidi->Fill(pos[0],pos[1]);
	255	uni->Fill(pos[2]);
	256	if(verbose){
	257	cout<<"==========================================================="<<endl;
	258	cout<<detna<<" module="<<mod<<endl;
	259	cout<<"Mod. coordinates: "<<pos[0]<<", "<<pos[1]<<", "<<pos[2]<<" Radius "<<ragdet<<endl;
	260	}
	261
	262	// Hits
	263	GetHitsCoor(ITS,mod,histos,subd,verbose);
	264
	265	//RecPoints
	266	if(userec){
	267	ITS->ResetRecPoints();
	268	TR->GetEvent(mod);
	269	TH2F bidi=(TH2F)histos.At(6+subd*9);
	270	TH1F uni=(TH1F)histos.At(7+subd*9);
	271	nrecp=GetRecCoor(geom,ITSrec,mod,bidi,uni,verbose);
	272	}
	273
	274	// Digits
	275	if(usedigits){
	276	ITS->ResetDigits();
	277	nbytes += TD->GetEvent(mod);
	278	GetDigits(seg,geom,ITSdigits,subd,mod,verbose,histos);
	279	}
	280
	281	} // End of loop on the modules
	282	TH1F *h1tmp;
	283	TH2F *h2tmp;
	284	// Plot the histograms
	285	TCanvas *current=0; // current Canvas (1--> SPD, 2---> SDD, 3---> SSD)
	286	if(subd==0){
	287	// Prepare canvas 1
	288	TCanvas *c1 = new TCanvas("c1","SPD",10,10,600,900);
	289	c1->Divide(2,3);
	290	current=c1;
	291	}
	292	if(subd==1){
	293	// Prepare canvas 2
	294	TCanvas *c2 = new TCanvas("c2","SDD",40,40,600,900);
	295	c2->Divide(2,3);
	296	current=c2;
	297	}
298	if(subd==2){
299	// Prepare canvas 3
300	TCanvas *c3 = new TCanvas("c3","SSD",70,70,600,900);
301	c3->Divide(2,3);
302	current=c3;
303	}
304	current->cd(1);
305	h2tmp = (TH2F)histos.At(9subd);
306	h2tmp->Draw();
307	current->cd(2);
308	h1tmp=(TH1F)histos.At(3+subd9);
309	h1tmp->Draw();
310	current->cd(3);
311	h2tmp=(TH2F)histos.At(1+9subd);
312	h2tmp->Draw();
313	current->cd(4);
314	h1tmp=(TH1F)histos.At(4+subd9);
315	h1tmp->Draw();
316
317	if(userec){
318	current->cd(5);
319	h2tmp=(TH2F)histos.At(6+9subd);
320	h2tmp->Draw();
321	current->cd(6);
322	h1tmp=(TH1F)histos.At(7+subd9);
323	h1tmp->Draw();
324	}
325
326	else {
327	current->cd(5);
328	h2tmp=(TH2F)histos.At(2+9subd);
329	h2tmp->Draw();
330	current->cd(6);
331	h2tmp=(TH2F)histos.At(5+9subd);
332	h2tmp->Draw();
333	}
334	} // if(All.....
335	} // end of loop on subdetectors
336	// Save the histograms
337	TFile *fh = new TFile("ITSgeoplot.root","recreate");
338	// The list is written to file as a single entry
339	TList *lihi = new TList();
340	// copy the pointers to the histograms to a TList object.
341	// The histograms concerning recpoints are not copied if
342	// 'userec' is false.
343	for(Int_t i=0;i<histos.GetEntriesFast();i++){
344	if(choice.Contains("All") \|\| (choice.Contains("SPD") && i<8) \|\| (choice.Contains("SDD") && i>7 && i<16) \|\| (choice.Contains("SSD") && i>15)){
345	if(!(!userec && ((i+2)%9==0 \|\| (i+1)%9==0)))lihi->Add(histos.At(i));
346	}
347	}
348	lihi->Write("Histograms ITS hits+digits+recpoints",TObject::kSingleKey);
349	fh->Close();
350
351	return retcode;
352	}
353
354
355	void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb){
356	TH2F h2=(TH2F)histos.At(1+subd*9);
357	TH1F h1=(TH1F)histos.At(4+subd*9);
358	TH1F ener=(TH1F)histos.At(8+subd*9);
359	AliITS ITS= (AliITS)its;
360	AliITSmodule *modu = ITS->GetModule(mod);
361	TObjArray *fHits = modu->GetHits();
362	Int_t nhits = fHits->GetEntriesFast();
363	if(nhits>0){
364	if(verb){
365	cout<<"-------------------------------------------------------"<<endl;
366	cout<<"Number of HITS for module "<<mod<<": "<<nhits<<endl;
367	}
368	for (Int_t hit=0;hit<nhits;hit++) {
369	AliITShit iHit = (AliITShit) fHits->At(hit);
370	if(!iHit->StatusEntering()){
371	Float_t x=iHit->GetXG();
372	Float_t y=iHit->GetYG();
373	Float_t z=iHit->GetZG();
374	Float_t edep=iHit->GetIonization()*1000000;
375	h2->Fill(x,y);
376	h1->Fill(z);
377	ener->Fill(edep);
378	if(verb){
379	cout<<"hit # "<<hit<<" Global coordinates "<<x<<" "<<y<<" "<<z<<endl;
380	cout<<"track # "<<iHit->GetTrack()<<" energy deposition (KeV)= "<<edep<<endl;
381	}
382	}
383	}
384	}
385	}
386
387
388	Int_t GetRecCoor(TObject ge, TClonesArray ITSrec, Int_t mod, TH2F h2, TH1F h1, Bool_t verb){
389	AliITSgeom geom = (AliITSgeom)ge;
390	Int_t nrecp = ITSrec->GetEntries();
391	if(nrecp>0){
392	Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
393	Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
394	if(verb){
395	cout<<"-------------------------------------------------------"<<endl;
396	cout<<"Number of REC POINTS for module "<<mod<<": "<<nrecp<<endl;
397	}
398	for(Int_t irec=0;irec<nrecp;irec++) {
399	AliITSRecPoint recp = (AliITSRecPoint)ITSrec->At(irec);
400	lc[0]=recp->GetX();
401	lc[2]=recp->GetZ();
402	geom->LtoG(mod,lc,gc);
403	if(verb){
404	cout<<"recp # "<<irec<<" local coordinates. lx= "<<lc[0]<<" lz= "<<lc[2]<<endl;
405	Float_t r=sqrt(gc[0]gc[0]+gc[1]gc[1]);
406	cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<r<<endl;
407	}
408	h2->Fill(gc[0],gc[1]);
409	h1->Fill(gc[2]);
410	}
411	}
412	return nrecp;
413	}
414
415	void GetDigits(TObject tmps,TObject ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos){
416	AliITSsegmentation seg = (AliITSsegmentation)tmps;
417	AliITSgeom geom = (AliITSgeom)ge;
79f97fd2	418	Int_t layer;
	419	Int_t ladder;
	420	Int_t detec;
	421	if(subd==2){
	422	geom->GetModuleId(mod,layer,ladder,detec);
	423	seg->SetLayer(layer);
	424	}
2f101c6a	425	Float_t lcoor[3]; for(Int_t j=0; j<3; j++) lcoor[j]=0.; //local coord dig.
	426	Float_t gcoor[3]; for(Int_t j=0; j<3; j++) gcoor[j]=0.; // global coo. dig.
	427	Float_t ragdig; // Radius digit
	428	TArrayI ssdone(5000); // used to match p and n side digits of strips
	429	TArrayI pair(5000); // as above
	430	Int_t ndigits = ITSdigits->GetEntries();
	431	AliITSdigit *digs;
	432	if(ndigits){
	433	if(verbose){
	434	cout<<"-------------------------------------------------------"<<endl;
	435	cout<<"Number of DIGITS for module "<<mod<<": "<<ndigits<<endl;
	436	}
	437	// Get the coordinates of the module
	438	if(subd==2){
	439	for (Int_t digit=0;digit<ndigits;digit++){
	440	ssdone[digit]=0;
	441	pair[digit]=0;
	442	}
	443	}
	444	for (Int_t digit=0;digit<ndigits;digit++) {
	445	digs = (AliITSdigit*)ITSdigits->UncheckedAt(digit);
	446	Int_t iz=digs->fCoord1; // cell number z
	447	Int_t ix=digs->fCoord2; // cell number x
faafc0eb	448	// Get local coordinates of the element
2f101c6a	449	if(subd<2){
faafc0eb	450	seg->DetToLocal(ix,iz,lcoor[0],lcoor[2]);
2f101c6a	451	}
	452	else{
	453	// SSD: if iz==0 ---> N side; if iz==1 P side
	454	if(ssdone[digit]==0){
	455	ssdone[digit]=1;
	456	pair[digit]=-1;
	457	Bool_t pside=(iz==1);
	458	Bool_t impaired=kTRUE;
	459	Int_t pstrip=0;
	460	Int_t nstrip=0;
	461	if(pside)pstrip=ix;
	462	if(!pside)nstrip=ix;
	463	for(Int_t digi2=0;digi2<ndigits;digi2++){
	464	if(ssdone[digi2]==0 && impaired){
	465	AliITSdigitSSD dig2=(AliITSdigitSSD)ITSdigits->UncheckedAt(digi2);
3a48961e	466	if(dig2->fCoord1 != iz && dig2->GetTrack(0)==digs->GetTrack(0) && dig2->GetTrack(0)>0){
2f101c6a	467	ssdone[digi2]=2;
	468	pair[digit]=digi2;
	469	if(pside)nstrip=dig2->fCoord2;
	470	if(!pside)pstrip=dig2->fCoord2;
	471	impaired=kFALSE;
	472	}
	473	}
	474	}
	475	if(!impaired)seg->GetPadCxz(pstrip,nstrip,lcoor[0],lcoor[2]);
	476	}
	477	}
2f101c6a	478	if(subd<2 \|\| (subd==2 && ssdone[digit]==1)){
	479	Int_t coor1=digs->fCoord1;
	480	Int_t coor2=digs->fCoord2;
3a48961e	481	Int_t tra0=digs->GetTrack(0);
2f101c6a	482	if(verbose){
	483	cout<<"digit # "<<digit<<" fCoord1= "<<coor1<<" fCoord2= "<<coor2<<" track "<<tra0<<" "<<endl;
	484	if(subd<2)cout<<"local coordinates -- x="<<lcoor[0]<<", z="<<lcoor[2]<<endl;
	485	if(subd==2){
	486	if(pair[digit]==-1){
	487	cout<<"(digit not paired)"<<endl;
	488	}
	489	else {
	490	cout<<"local coordinates -- x="<<lcoor[0]<<", z="<<lcoor[2]<<endl;
	491	Int_t dtmp=pair[digit];
	492	AliITSdigitSSD dig2=(AliITSdigitSSD)ITSdigits->UncheckedAt(dtmp);
	493	Int_t coor1b=dig2->fCoord1;
	494	Int_t coor2b=dig2->fCoord2;
3a48961e	495	Int_t tra0b=dig2->GetTrack(0);
2f101c6a	496	cout<<"(digit paired with digit #"<<dtmp<<endl;
	497	cout<<"with fCoord1= "<<coor1b<<" fCoord2= "<<coor2b<<" track "<<tra0b<<")"<<endl;
	498	}
	499	}
	500	}
	501	if(subd<2 \|\| (subd==2 && pair[digit]!=-1)){
	502	// Global coordinates of the element
faafc0eb	503	//SDD and SPD use cm, SSD microns (GetPadCxz)
faafc0eb	504	if(subd==2)for(Int_t j=0;j<3;j++)lcoor[j]=lcoor[j]/10000.;
2f101c6a	505	lcoor[1]=0.;
	506	geom->LtoG(mod,lcoor,gcoor); // global coord. in cm
	507	ragdig=sqrt(gcoor[0]gcoor[0]+gcoor[1]gcoor[1]);
	508	if(verbose)cout<<"global coordinates "<<gcoor[0]<<" "<<gcoor[1]<<" "<<gcoor[2]<<" Radius "<<ragdig<<endl;
	509	//Fill histograms
	510	TH2F bidi = (TH2F)histos.At(subd*9);
	511	TH1F uni = (TH1F)histos.At(3+subd*9);
	512	bidi->Fill(gcoor[0],gcoor[1]);
	513	uni->Fill(gcoor[2]);
	514	}
	515	}
	516	} // loop on digits for this module
	517	} // if(ndigits>0....
	518	}
	519
	520
	521
	522
	523
	524
	525
	526
	527
	528
	529
	530
	531
	532
	533
	534
	535
	536
	537
	538