ITS/AliITSGeoPlot.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include<Riostream.h>
   3 #include<TROOT.h>
   4 #include<TArrayI.h>
   5 #include<TBranch.h>
   6 #include<TCanvas.h>
   7 #include<TClassTable.h>
   8 #include<TClonesArray.h>
   9 #include<TFile.h>
  10 #include<TGeoManager.h>
  11 #include<TH1.h>
  12 #include<TH2.h>
  13 #include <TInterpreter.h>
  14 #include<TObject.h>
  15 #include<TObjArray.h>
  16 #include<TTree.h>
  17 #include "AliRun.h"
  18 #include "AliITS.h"
  19 #include "AliITSgeom.h"
  20 #include "AliITSDetTypeRec.h"
  21 #include "AliITSRecPoint.h"
  22 #include "AliITSRecPoint.h"
  23 #include "AliITSdigit.h"
  24 #include "AliITSdigitSSD.h"
  25 #include "AliITShit.h"
  26 #include "AliITSmodule.h"
  27 #include "AliITSsegmentation.h"
  28 #include "AliITSsegmentationSPD.h"
  29 #include "AliITSsegmentationSDD.h"
  30 #include "AliITSsegmentationSSD.h"
  31 #include "AliRunLoader.h"
  32 #include "AliITSLoader.h"
  33 #include "AliHeader.h"
  34 #include "AliCDBManager.h"
  35 #include "AliCDBStorage.h"
  36 #endif
  37 void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb);
  38 Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
  39 Int_t GetClusCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
  40 void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos);
  41
  42 Int_t AliITSGeoPlot (Int_t evesel=0, char *opt="All+ClustersV2", TString filename="galice.root", Int_t isfastpoints = 0) {
  43   /*******************************************************************
  44    *  This macro displays geometrical information related to the
  45    *  hits, digits and rec points (or V2 clusters) in ITS.
  46    *  There are histograms that are not displayed (i.e. energy
  47    *  deposition) but that are saved onto a file (see below)
  48    *
  49    *  INPUT arguments:
  50    *
  51    *  Int_t evesel:  the selected event number. It's =0 by default
  52    *
  53    *  Options: Any combination of:
  54    *    1) subdetector name:  "SPD", "SDD", "SSD", "All" (default)
  55    *    2) Printouts:        "Verbose" Almost everything is printed out:
  56    *                          it is wise to redirect the output onto a file
  57    *                    e.g.: .x AliITSGeoPlot.C("All+Verbose") > out.log
  58    *    3) Rec Points option: "Rec"   ---> Uses Rec. Points (default)
  59    *
  60    *    4) ClustersV2 option: "ClustersV2" ---->Uses ClustersV2
  61    *                           otherwise ---->uses hits and digits only
  62    *    Examples:
  63    *       .x AliITSGeoPlot();  (All subdetectors; no-verbose; no-recpoints)
  64    *       .x AliITSGeoPlot("SPD+SSD+Verbose+Rec");
  65    *
  66    *    filename:   It's "galice.root" by default.
  67    *    isfastpoints: integer. It is set to 0 by defaults. This means that
  68    *                slow reconstruction is assumed. If fast recpoint are in
  69    *                in use, isfastpoints must be set =1.
  70    *
  71    *  OUTPUT: It produces a root file with a list of histograms
  72    *
  73    *  WARNING: spatial information for SSD/DIGITS is obtained by pairing
  74    *          digits on p and n side originating from the same track, when
  75    *          possible. This (mis)use of DIGITS is tolerated for debugging
  76    *          purposes only !!!!  The pairing in real life should be done
  77    *          starting from info really available...
  78    *
  79    *  COMPILATION: this macro can be compiled.
  80    *      1)       You need to set your include path with
  81    * gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g");
  82    *      3)       If you are using root instead of aliroot you need to
  83    *               execute the macro loadlibs.C in advance
  84    *      4)       To compile this macro from root (or aliroot):
  85    *                 ---  .L AliITSGeoPlot.C++
  86    *                 ---  AliITSGeoPlot();
  87    *
  88    *  M.Masera  14/05/2001 18:30
  89    *  Last rev. 31/05/2004 14:00 (Clusters V2 added)  E.C.
  90    ********************************************************************/
  91
  92   //Options
  93   TString choice(opt);
  94   Bool_t All = choice.Contains("All");
  95   Bool_t verbose=choice.Contains("Verbose");
  96   Bool_t userec=choice.Contains("Rec");
  97   Bool_t useclustersv2=choice.Contains("ClustersV2");
  98   Int_t retcode=1; //return code
  99   if (gClassTable->GetID("AliRun") < 0) {
 100     gInterpreter->ExecuteMacro("loadlibs.C");
 101   }
 102   else {
 103     if(gAlice){
 104       delete gAlice->GetRunLoader();
 105       delete gAlice;
 106       gAlice=0;
 107     }
 108   }
 109   // Set OCDB if needed
 110   AliCDBManager* man = AliCDBManager::Instance();
 111   if (!man->IsDefaultStorageSet()) {
 112     printf("Setting a local default storage and run number 0\n");
 113     man->SetDefaultStorage("local://$ALICE_ROOT");
 114     man->SetRun(0);
 115   }
 116   else {
 117     printf("Using deafult storage \n");
 118   }
 119   // retrives geometry
 120   TString geof(gSystem->DirName(filename));
 121   geof += "/geometry.root";
 122   TGeoManager::Import(geof.Data());
 123   if (!gGeoManager) {
 124     cout<<"geometry not found\n";
 125     return -1;
 126   }
 127
 128   AliRunLoader* rl = AliRunLoader::Open(filename.Data());
 129   if (rl == 0x0){
 130     cerr<<"AliITSGeoPlot.C : Can not open session RL=NULL"<< endl;
 131     return -1;
 132   }
 133   Int_t retval = rl->LoadgAlice();
 134   if (retval){
 135     cerr<<"AliITSGeoPlot.C : LoadgAlice returned error"<<endl;
 136     return -1;
 137   }
 138   gAlice=rl->GetAliRun();
 139
 140   retval = rl->LoadHeader();
 141   if (retval){
 142     cerr<<"AliITSGeoPlot.C : LoadHeader returned error"<<endl;
 143     return -1;
 144   }
 145
 146   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 147
 148   if(!ITSloader){
 149     cerr<<"AliITSGeoPlot.C :  ITS loader not found"<<endl;
 150     return -1;
 151   }
 152
 153   ITSloader->LoadHits("read");
 154   ITSloader->LoadDigits("read");
 155   if(isfastpoints==1)ITSloader->SetRecPointsFileName("ITS.FastRecPoints.root");
 156   ITSloader->LoadRecPoints("read");
 157   rl->GetEvent(evesel);
 158   Int_t nparticles = rl->GetHeader()->GetNtrack();
 159   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 160   ITS->SetTreeAddress();
 161   if(verbose) {
 162     cout<<" "<<endl<<" "<<endl;
 163     cout<<"******* Event processing started   *******"<<endl;
 164     cout<<"In the following, hits with 'StatusEntering' flag active"<<endl;
 165     cout<<"will not be processed"<<endl;
 166     cout << "Number of particles=  " << nparticles <<endl;
 167   }
 168
 169   // HITS
 170   TTree *TH = ITSloader->TreeH();
 171   Stat_t ntracks = TH->GetEntries();
 172   if(verbose)cout<<"Number of primary tracks= "<<ntracks<<endl;
 173
 174   // ITS
 175   Int_t nmodules;
 176   ITS->InitModules(-1,nmodules);
 177   cout<<"Number of ITS modules= "<<nmodules<<endl;
 178   cout<<"Filling modules... It takes a while, now. Please be patient"<<endl;
 179   ITS->FillModules(0,0,nmodules," "," ");
 180   cout<<"ITS modules .... DONE!"<<endl;
 181
 182   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec(ITSloader);
 183   detTypeRec->SetDefaults();
 184   // DIGITS
 185   TTree *TD = ITSloader->TreeD();
 186
 187   //RECPOINTS (V2 clusters)
 188   TTree *TR = ITSloader->TreeR();
 189   TClonesArray *ITSrec  = detTypeRec->RecPoints();
 190   TBranch *branch = 0;
 191   if(userec && TR && ITSrec){
 192     if(isfastpoints==1){
 193       branch = ITSloader->TreeR()->GetBranch("ITSRecPointsF");
 194       cout<<"using fast points\n";
 195     }
 196     else {
 197       branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
 198     }
 199     if(branch)branch->SetAddress(&ITSrec);
 200   }
 201
 202   if(userec && (!TR || !ITSrec || !branch)){
 203     userec = kFALSE;
 204     cout<<"\n ======================================================= \n";
 205     cout<<"WARNING: there are no RECPOINTS on this file ! \n";
 206     cout<<"======================================================= \n \n";
 207   }
 208   if(useclustersv2 && TR && ITSrec){
 209     branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
 210     if(branch)branch->SetAddress(&ITSrec);
 211   }
 212
 213   if(useclustersv2 && (!TR || !ITSrec || !branch)){
 214     useclustersv2 = kFALSE;
 215     cout<<"\n ======================================================= \n";
 216     cout<<"WARNING: there are no CLUSTERSV2 on this file ! \n";
 217     cout<<"======================================================= \n \n";
 218   }
 219
 220
 221   //local variables
 222   Int_t mod;   //module number
 223   Int_t nbytes = 0;
 224   Double_t pos[3];  // Global position of the current module
 225   Float_t ragdet; // Radius of detector (x y plane)
 226   Int_t first; // first module
 227   Int_t last; // last module
 228   Int_t nrecp; //number of RecPoints for a given module
 229
 230   //List of histograms
 231   TObjArray histos(26,0);  // contains the pointers to the histograms
 232   // Book histograms SPD
 233   TH2F *rspd = new TH2F("rspd","Radii of digits - SPD",50,-10.,10.,50,-10.,10.);
 234   TH2F *rhspd = new TH2F("rhspd","Radii of hits - SPD",50,-10.,10.,50,-10.,10.);
 235   TH2F *rmspd = new TH2F("rmspd","Radii of SPD modules",50,-10.,10.,50,-10.,10.);
 236   TH1F *zspd = new TH1F("zspd","Z of digits - SPD",100,-30.,30.);
 237   TH1F *zhspd = new TH1F("zhspd","Z of hits - SPD",100,-30.,30.);
 238   TH1F *zmspd = new TH1F("zmspd","Z of SPD modules",100,-30,30.);
 239
 240   Char_t title1[50]="";
 241   Char_t title2[50]="";
 242   if(userec){
 243     sprintf(title1,"Radii of recpoints - %s","SPD");
 244     sprintf(title2,"Z of recpoints - %s","SPD");
 245   }
 246   if(useclustersv2){
 247     sprintf(title1,"Radii of clustersV2 - %s","SPD");
 248     sprintf(title2,"Z of clustersV2 - %s","SPD");
 249   }
 250   TH2F *rrspd = new TH2F("rrspd",title1,50,-10.,10.,50,-10.,10.);
 251   TH1F *zrspd = new TH1F("zrspd",title2,100,-30.,30.);
 252   TH1F *enespd = new TH1F("enespd","Energy deposition SPD (KeV)",100,0.,1000.);
 253   histos.AddLast(rspd);  // 0
 254   histos.AddLast(rhspd); // 1
 255   histos.AddLast(rmspd); // 2
 256   histos.AddLast(zspd);  // 3
 257   histos.AddLast(zhspd); // 4
 258   histos.AddLast(zmspd); // 5
 259   histos.AddLast(rrspd); // 6
 260   histos.AddLast(zrspd); // 7
 261   histos.AddLast(enespd); // 8
 262   // Book histograms SDD
 263   TH2F *rsdd = new TH2F("rsdd","Radii of digits - SDD",50,-40.,40.,50,-40.,40.);
 264   TH2F *rhsdd = new TH2F("rhsdd","Radii of hits - SDD",50,-40.,40.,50,-40.,40.);
 265   TH2F *rmsdd = new TH2F("rmsdd","Radii of SDD modules",50,-40.,40.,50,-40.,40.);
 266   TH1F *zsdd = new TH1F("zsdd","Z of digits - SDD",100,-40.,40.);
 267   TH1F *zhsdd = new TH1F("zhsdd","Z of hits - SDD",100,-40.,40.);
 268   TH1F *zmsdd = new TH1F("zmsdd","Z of SDD modules",100,-40,40.);
 269   Char_t title3[50];
 270   Char_t title4[50];
 271   if(userec){
 272     sprintf(title3,"Radii of recpoints - %s","SDD");
 273     sprintf(title4,"Z of recpoints - %s","SDD");
 274   }
 275   if(useclustersv2){
 276     sprintf(title3,"Radii of clustersV2 - %s","SDD");
 277     sprintf(title4,"Z of clustersV2 - %s","SDD");
 278   }
 279   TH2F *rrsdd = new TH2F("rrsdd",title3,50,-40.,40.,50,-40.,40.);
 280   TH1F *zrsdd = new TH1F("zrsdd",title4,100,-40.,40.);
 281   TH1F *enesdd = new TH1F("enesdd","Energy deposition SDD (KeV)",100,0.,1000.);
 282   histos.AddLast(rsdd);  // 9
 283   histos.AddLast(rhsdd); // 10
 284   histos.AddLast(rmsdd); // 11
 285   histos.AddLast(zsdd);  // 12
 286   histos.AddLast(zhsdd); // 13
 287   histos.AddLast(zmsdd); // 14
 288   histos.AddLast(rrsdd); // 15
 289   histos.AddLast(zrsdd); // 16
 290   histos.AddLast(enesdd); // 17
 291   // Book histogram SSD
 292   TH2F *rssd = new TH2F("rssd","Radii of digits - SSD",50,-50.,50.,50,-50.,50.);
 293   TH2F *rhssd = new TH2F("rhssd","Radii of hits - SSD",50,-50.,50.,50,-50.,50.);
 294   TH2F *rmssd = new TH2F("rmssd","Radii of SSD modules",50,-50.,50.,50,-50.,50.);
 295   TH1F *zssd = new TH1F("zssd","Z of digits - SSD",100,-70.,70.);
 296   TH1F *zhssd = new TH1F("zhssd","Z of hits - SSD",100,-70.,70.);
 297   TH1F *zmssd = new TH1F("zmssd","Z of SSD modules",100,-70,70.);
 298   Char_t title5[50];
 299   Char_t title6[50];
 300   if(userec){
 301     sprintf(title5,"Radii of recpoints - %s","SSD");
 302     sprintf(title6,"Z of recpoints - %s","SSD");
 303   }
 304   if(useclustersv2){
 305     sprintf(title5,"Radii of clustersV2 - %s","SSD");
 306     sprintf(title6,"Z of clustersV2 - %s","SSD");
 307   }
 308
 309   TH2F *rrssd = new TH2F("rrssd",title5,50,-50.,50.,50,-50.,50.);
 310   TH1F *zrssd = new TH1F("zrssd",title6,100,-70.,70.);
 311   TH1F *enessd = new TH1F("enessd","Energy deposition SSD (KeV)",100,0.,1000.);
 312   histos.AddLast(rssd);  // 18
 313   histos.AddLast(rhssd); // 19
 314   histos.AddLast(rmssd); // 20
 315   histos.AddLast(zssd);  // 21
 316   histos.AddLast(zhssd); // 22
 317   histos.AddLast(zmssd); // 23
 318   histos.AddLast(rrssd); // 24
 319   histos.AddLast(zrssd); // 25
 320   histos.AddLast(enessd); // 26
 321   //
 322   // Loop on subdetectors
 323   //
 324   cout<<"CALL GETITSGEOM \n \n \n";
 325   AliITSgeom *geom = ITS->GetITSgeom();
 326   TString detna; // subdetector name
 327   for(Int_t subd=0;subd<3;subd++){
 328     if(All || (choice.Contains("SPD") && subd==0) || (choice.Contains("SDD") && subd==1) || (choice.Contains("SSD") && subd==2)){
 329       // Prepare array for the digits
 330       TClonesArray *ITSdigits  = ITS->DigitsAddress(subd);
 331       Bool_t usedigits = kTRUE;
 332       if(!ITSdigits){
 333         usedigits = kFALSE;
 334         cout<<"\n ======================================================= \n";
 335         cout<<"WARNING: there are no DIGITS on this file ! \n";
 336         cout<<"======================================================= \n \n";
 337       }
 338       // Get segmentation model
 339       if(subd==0)detna="SPD";
 340       if(subd==1)detna="SDD";
 341       if(subd==2)detna="SSD";
 342       AliITSsegmentation *seg=(AliITSsegmentation*)detTypeRec->GetSegmentationModel(subd);
 343       // Loop on modules
 344       first = geom->GetStartDet(subd);
 345       last = geom->GetLastDet(subd);
 346       if(verbose){
 347         cout<<"     "<<endl<<"-------------------------------------"<<endl;
 348         cout<<"Start processing subdetector "<<detna<<endl;
 349         cout<<detna<<" first module "<<first<<endl;
 350         cout<<detna<<" last module "<<last<<endl;
 351         cout<<" "<<endl<<" "<<endl;
 352       }
 353       for (mod=first; mod<=last; mod++){
 354         geom->GetTrans(mod,pos);  // position of the module in the MRS
 355         ragdet=sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
 356         // The following 2 histos are a check of the geometry
 357         TH2F *bidi = (TH2F*)histos.At(2+subd*9);
 358         TH1F *uni = (TH1F*)histos.At(5+subd*9);
 359         bidi->Fill(pos[0],pos[1]);
 360         uni->Fill(pos[2]);
 361         if(verbose){
 362           cout<<"=========================================================\n";
 363           cout<<detna<<" module="<<mod<<endl;
 364           cout<<"Mod. coordinates: "<<pos[0]<<", "<<pos[1]<<", ";
 365           cout<<pos[2]<<" Radius "<<ragdet<<endl;
 366         }
 367
 368         // Hits
 369         GetHitsCoor(ITS,mod,histos,subd,verbose);
 370
 371         //RecPoints
 372         if(userec){
 373           detTypeRec->ResetRecPoints();
 374           branch->GetEvent(mod);
 375           TH2F *bidi=(TH2F*)histos.At(6+subd*9);
 376           TH1F *uni=(TH1F*)histos.At(7+subd*9);
 377           nrecp=GetRecCoor(geom,ITSrec,mod,bidi,uni,verbose);
 378         }
 379         if(useclustersv2){
 380           detTypeRec->ResetRecPoints();
 381           branch->GetEvent(mod);
 382           TH2F *bidi=(TH2F*)histos.At(6+subd*9);
 383           TH1F *uni=(TH1F*)histos.At(7+subd*9);
 384           nrecp=GetClusCoor(geom,ITSrec,mod,bidi,uni,verbose);
 385
 386         }
 387
 388         // Digits
 389         if(usedigits){
 390           detTypeRec->ResetDigits();
 391           nbytes += TD->GetEvent(mod);
 392           GetDigits(seg,geom,ITSdigits,subd,mod,verbose,histos);
 393         }
 394
 395       } // End of loop on the modules
 396       TH1F *h1tmp;
 397       TH2F *h2tmp;
 398       //  Plot the histograms
 399       TCanvas *current=0; // current Canvas (1--> SPD, 2---> SDD, 3---> SSD)
 400       if(subd==0){
 401         // Prepare canvas 1
 402         TCanvas *c1 = new TCanvas("c1","SPD",10,10,600,900);
 403         c1->Divide(2,3);
 404         current=c1;
 405       }
 406       if(subd==1){
 407         // Prepare canvas 2
 408         TCanvas *c2 = new TCanvas("c2","SDD",40,40,600,900);
 409         c2->Divide(2,3);
 410         current=c2;
 411       }
 412       if(subd==2){
 413         // Prepare canvas 3
 414         TCanvas *c3 = new TCanvas("c3","SSD",70,70,600,900);
 415         c3->Divide(2,3);
 416         current=c3;
 417       }
 418       current->cd(1);
 419       h2tmp = (TH2F*)histos.At(9*subd);
 420       h2tmp->Draw();
 421       current->cd(2);
 422       h1tmp=(TH1F*)histos.At(3+subd*9);
 423       h1tmp->Draw();
 424       current->cd(3);
 425       h2tmp=(TH2F*)histos.At(1+9*subd);
 426       h2tmp->Draw();
 427       current->cd(4);
 428       h1tmp=(TH1F*)histos.At(4+subd*9);
 429       h1tmp->Draw();
 430
 431       if(userec || useclustersv2){
 432         current->cd(5);
 433         h2tmp=(TH2F*)histos.At(6+9*subd);
 434         h2tmp->Draw();
 435         current->cd(6);
 436         h1tmp=(TH1F*)histos.At(7+subd*9);
 437         h1tmp->Draw();
 438       }
 439
 440       else {
 441         current->cd(5);
 442         h2tmp=(TH2F*)histos.At(2+9*subd);
 443         h2tmp->Draw();
 444         current->cd(6);
 445         h2tmp=(TH2F*)histos.At(5+9*subd);
 446         h2tmp->Draw();
 447       }
 448     } // if(All.....
 449   } // end of loop on subdetectors
 450   // Save the histograms
 451   TFile *fh = new TFile("AliITSGeoPlot.root","recreate");
 452   // The list is written to file as a single entry
 453   TList *lihi = new TList();
 454   // copy the pointers to the histograms to a TList object.
 455   // The histograms concerning recpoints are not copied if
 456   // 'userec' is false.
 457   for(Int_t i=0;i<histos.GetEntriesFast();i++){
 458     if(choice.Contains("All") || (choice.Contains("SPD") && i<8) || (choice.Contains("SDD") && i>7 && i<16) || (choice.Contains("SSD") && i>15)){
 459       if(!(!userec && ((i+2)%9==0 || (i+1)%9==0)))lihi->Add(histos.At(i));
 460     }
 461   }
 462   lihi->Write("Histograms ITS hits+digits+recpoints",TObject::kSingleKey);
 463   fh->Close();
 464
 465   return retcode;
 466 }
 467
 468
 469 void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb){
 470   TH2F *h2=(TH2F*)histos.At(1+subd*9);
 471   TH1F *h1=(TH1F*)histos.At(4+subd*9);
 472   TH1F *ener=(TH1F*)histos.At(8+subd*9);
 473   AliITS *ITS= (AliITS*)its;
 474   AliITSmodule *modu = ITS->GetModule(mod);
 475   TObjArray *fHits = modu->GetHits();
 476   Int_t nhits = fHits->GetEntriesFast();
 477   if(nhits>0){
 478     if(verb){
 479       cout<<"-------------------------------------------------------"<<endl;
 480       cout<<"Number of HITS for module "<<mod<<": "<<nhits<<endl;
 481     }
 482     for (Int_t hit=0;hit<nhits;hit++) {
 483       AliITShit *iHit = (AliITShit*) fHits->At(hit);
 484       if(!iHit->StatusEntering()){
 485         Float_t x=iHit->GetXG();
 486         Float_t y=iHit->GetYG();
 487         Float_t z=iHit->GetZG();
 488         Float_t edep=iHit->GetIonization()*1000000;
 489         h2->Fill(x,y);
 490         h1->Fill(z);
 491         ener->Fill(edep);
 492         if(verb){
 493           cout<<"hit # "<<hit<<" Global coordinates "<<x<<" "<<y<<" "<<z<<endl;
 494           cout<<"track # "<<iHit->GetTrack()<<" energy deposition (KeV)= ";
 495           cout<<edep<<endl;
 496         }
 497       }
 498     }
 499   }
 500 }
 501
 502
 503 Int_t GetClusCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
 504
 505   AliITSgeom *geom = (AliITSgeom*)ge;
 506   Int_t nrecp = ITSrec->GetEntries();
 507   if(nrecp>0){
 508     Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
 509     Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
 510     if(verb){
 511       cout<<"-------------------------------------------------------"<<endl;
 512       cout<<"Number of CLUSTERS for module "<<mod<<": "<<nrecp<<endl;
 513     }
 514     for(Int_t irec=0;irec<nrecp;irec++) {
 515       AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
 516       Double_t rot[9];
 517       geom->GetRotMatrix(mod,rot);
 518       Int_t lay,lad,det;
 519       geom->GetModuleId(mod,lay,lad,det);
 520       Float_t tx,ty,tz;
 521       geom->GetTrans(lay,lad,det,tx,ty,tz);
 522
 523       Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
 524       Double_t phi1=TMath::Pi()/2+alpha;
 525       if(lay==1) phi1+=TMath::Pi();
 526
 527       Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
 528       Float_t  r=tx*cp+ty*sp;
 529       gc[0]= r*cp - recp->GetY()*sp;
 530       gc[1]= r*sp + recp->GetY()*cp;
 531       gc[2]=recp->GetZ();
 532
 533       if(verb){
 534         Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
 535         cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" ";
 536         cout<<r<<endl;
 537         cout<<"Associated track "<<recp->GetLabel(0)<<endl;
 538       }
 539       h2->Fill(gc[0],gc[1]);
 540       h1->Fill(gc[2]);
 541
 542     }
 543   }
 544   return nrecp;
 545 }
 546 Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
 547   AliITSgeom *geom = (AliITSgeom*)ge;
 548   Int_t nrecp = ITSrec->GetEntries();
 549   if(nrecp>0){
 550     Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
 551     Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
 552     if(verb){
 553       cout<<"-------------------------------------------------------"<<endl;
 554       cout<<"Number of REC POINTS for module "<<mod<<": "<<nrecp<<endl;
 555     }
 556     for(Int_t irec=0;irec<nrecp;irec++) {
 557       AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
 558       lc[0]=recp->GetDetLocalX();
 559       lc[2]=recp->GetDetLocalZ();
 560       geom->LtoG(mod,lc,gc);
 561       if(verb){
 562         cout<<"recp # "<<irec<<" local coordinates. lx= "<<lc[0]<<" lz= ";
 563         cout<<lc[2]<<endl;
 564         Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
 565         cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" ";
 566         cout<<r<<endl;
 567         cout<<"Associated track "<<recp->GetLabel(0)<<endl;
 568       }
 569       h2->Fill(gc[0],gc[1]);
 570       h1->Fill(gc[2]);
 571     }
 572   }
 573   return nrecp;
 574 }
 575
 576 void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos){
 577   AliITSsegmentation *seg = (AliITSsegmentation*)tmps;
 578   AliITSgeom *geom = (AliITSgeom*)ge;
 579   Int_t layer;
 580   Int_t ladder;
 581   Int_t detec;
 582   if(subd==2){
 583     geom->GetModuleId(mod,layer,ladder,detec);
 584     seg->SetLayer(layer);
 585   }
 586   Float_t lcoor[3]; for(Int_t j=0; j<3; j++) lcoor[j]=0.;  //local coord dig.
 587   Float_t gcoor[3]; for(Int_t j=0; j<3; j++) gcoor[j]=0.; // global coo. dig.
 588   Float_t ragdig; // Radius digit
 589   TArrayI ssdone(5000);  // used to match p and n side digits of strips
 590   TArrayI pair(5000);    // as above
 591   Int_t ndigits = ITSdigits->GetEntries();
 592   AliITSdigit *digs;
 593   if(ndigits){
 594     if(verbose){
 595       cout<<"-------------------------------------------------------"<<endl;
 596       cout<<"Number of DIGITS for module "<<mod<<": "<<ndigits<<endl;
 597     }
 598     // Get the coordinates of the module
 599     if(subd==2){
 600       for (Int_t digit=0;digit<ndigits;digit++){
 601         ssdone[digit]=0;
 602         pair[digit]=0;
 603       }
 604     }
 605     for (Int_t digit=0;digit<ndigits;digit++) {
 606       digs = (AliITSdigit*)ITSdigits->UncheckedAt(digit);
 607       Int_t iz=digs->GetCoord1();  // cell number z
 608       Int_t ix=digs->GetCoord2();  // cell number x
 609       // Get local coordinates of the element
 610       if(subd<2){
 611         seg->DetToLocal(ix,iz,lcoor[0],lcoor[2]);
 612       }
 613       else{
 614         // SSD: if iz==0 ---> N side; if iz==1 P side
 615         if(ssdone[digit]==0){
 616           ssdone[digit]=1;
 617           pair[digit]=-1;
 618           Bool_t pside=(iz==1);
 619           Bool_t impaired=kTRUE;
 620           Int_t pstrip=0;
 621           Int_t nstrip=0;
 622           if(pside)pstrip=ix;
 623           if(!pside)nstrip=ix;
 624           for(Int_t digi2=0;digi2<ndigits;digi2++){
 625             if(ssdone[digi2]==0 && impaired){
 626               AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(digi2);
 627               if(dig2->GetCoord1() != iz && dig2->GetTrack(0)==digs->GetTrack(0) && dig2->GetTrack(0)>0){
 628                 ssdone[digi2]=2;
 629                 pair[digit]=digi2;
 630                 if(pside)nstrip=dig2->GetCoord2();
 631                 if(!pside)pstrip=dig2->GetCoord2();
 632                 impaired=kFALSE;
 633               }
 634             }
 635           }
 636           if(!impaired)seg->GetPadCxz(pstrip,nstrip,lcoor[0],lcoor[2]);
 637         }
 638       }
 639       if(subd<2 || (subd==2 && ssdone[digit]==1)){
 640         Int_t coor1=digs->GetCoord1();
 641         Int_t coor2=digs->GetCoord2();
 642         Int_t tra0=digs->GetTrack(0);
 643         if(verbose){
 644           cout<<"digit # "<<digit<<" fCoord1= "<<coor1<<" fCoord2= ";
 645           cout<<coor2<<" track "<<tra0<<" "<<endl;
 646           if(subd<2)cout<<"local coordinates -- x="<<lcoor[0]<<", z=";
 647           cout<<lcoor[2]<<endl;
 648           if(subd==2){
 649             if(pair[digit]==-1){
 650               cout<<"(digit  not paired)"<<endl;
 651             }
 652             else {
 653               cout<<"local coordinates -- x="<<lcoor[0]<<", z="<<lcoor[2];
 654               cout<<endl;
 655               Int_t dtmp=pair[digit];
 656               AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(dtmp);
 657               Int_t coor1b=dig2->GetCoord1();
 658               Int_t coor2b=dig2->GetCoord2();
 659               Int_t tra0b=dig2->GetTrack(0);
 660               cout<<"(digit paired with digit #"<<dtmp<<endl;
 661               cout<<"with fCoord1= "<<coor1b<<" fCoord2= "<<coor2b;
 662               cout<<" track "<<tra0b<<")"<<endl;
 663             }
 664           }
 665         }
 666         if(subd<2 || (subd==2 && pair[digit]!=-1)){
 667           // Global coordinates of the element
 668           //SDD and SPD use cm, SSD microns (GetPadCxz)
 669           if(subd==2)for(Int_t j=0;j<3;j++)lcoor[j]=lcoor[j]/10000.;
 670           lcoor[1]=0.;
 671           geom->LtoG(mod,lcoor,gcoor);  // global coord. in cm
 672           ragdig=sqrt(gcoor[0]*gcoor[0]+gcoor[1]*gcoor[1]);
 673           if(verbose){
 674             cout<<"global coordinates "<<gcoor[0]<<" "<<gcoor[1];
 675             cout<<" "<<gcoor[2]<<" Radius "<<ragdig<<endl;
 676           }
 677           //Fill histograms
 678           TH2F *bidi = (TH2F*)histos.At(subd*9);
 679           TH1F *uni = (TH1F*)histos.At(3+subd*9);
 680           bidi->Fill(gcoor[0],gcoor[1]);
 681           uni->Fill(gcoor[2]);
 682         }
 683       }
 684     } // loop on digits for this module
 685   } // if(ndigits>0....
 686 }
 687