ITS/AliITSGeoPlot.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2
   3 #include<Riostream.h>
   4
   5 #include<TROOT.h>
   6
   7 #include<TArrayI.h>
   8
   9 #include<TBranch.h>
  10
  11 #include<TCanvas.h>
  12
  13 #include<TClassTable.h>
  14
  15 #include<TClonesArray.h>
  16
  17 #include<TFile.h>
  18
  19 #include<TH1.h>
  20
  21 #include<TH2.h>
  22
  23 #include <TInterpreter.h>
  24
  25 #include<TObject.h>
  26
  27 #include<TObjArray.h>
  28
  29 #include<TTree.h>
  30
  31 #include "AliRun.h"
  32
  33 #include "AliITS.h"
  34
  35 #include "AliITSgeom.h"
  36
  37 #include "AliITSDetTypeRec.h"
  38
  39 #include "AliITSRecPoint.h"
  40
  41 #include "AliITSRecPoint.h"
  42
  43 #include "AliITSdigit.h"
  44
  45 #include "AliITSdigitSSD.h"
  46
  47 #include "AliITShit.h"
  48
  49 #include "AliITSmodule.h"
  50
  51 #include "AliITSsegmentation.h"
  52
  53 #include "AliITSsegmentationSPD.h"
  54
  55 #include "AliITSsegmentationSDD.h"
  56
  57 #include "AliITSsegmentationSSD.h"
  58
  59 #include "AliRunLoader.h"
  60
  61 #include "AliITSLoader.h"
  62
  63 #include "AliHeader.h"
  64
  65 #endif
  66
  67 void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb);
  68
  69 Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
  70
  71 Int_t GetClusCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
  72
  73 void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos);
  74
  75
  76
  77 Int_t AliITSGeoPlot (Int_t evesel=0, char *opt="All+ClustersV2", char *filename="galice.root", Int_t isfastpoints = 0) {
  78
  79   /*******************************************************************
  80
  81    *  This macro displays geometrical information related to the
  82
  83    *  hits, digits and rec points (or V2 clusters) in ITS.
  84
  85    *  There are histograms that are not displayed (i.e. energy
  86
  87    *  deposition) but that are saved onto a file (see below)
  88
  89    *
  90
  91    *  INPUT arguments:
  92
  93    *
  94
  95    *  Int_t evesel:  the selected event number. It's =0 by default
  96
  97    *
  98
  99    *  Options: Any combination of:
 100
 101    *    1) subdetector name:  "SPD", "SDD", "SSD", "All" (default)
 102
 103    *    2) Printouts:        "Verbose" Almost everything is printed out:
 104
 105    *                          it is wise to redirect the output onto a file
 106
 107    *                    e.g.: .x AliITSGeoPlot.C("All+Verbose") > out.log
 108
 109    *    3) Rec Points option: "Rec"   ---> Uses Rec. Points (default)
 110
 111    *
 112
 113    *    4) ClustersV2 option: "ClustersV2" ---->Uses ClustersV2
 114
 115    *                           otherwise ---->uses hits and digits only
 116
 117    *    Examples:
 118
 119    *       .x AliITSGeoPlot();  (All subdetectors; no-verbose; no-recpoints)
 120
 121    *       .x AliITSGeoPlot("SPD+SSD+Verbose+Rec");
 122
 123    *
 124
 125    *    filename:   It's "galice.root" by default.
 126
 127    *    isfastpoints: integer. It is set to 0 by defaults. This means that
 128
 129    *                slow reconstruction is assumed. If fast recpoint are in
 130
 131    *                in use, isfastpoints must be set =1.
 132
 133    *
 134
 135    *  OUTPUT: It produces a root file with a list of histograms
 136
 137    *
 138
 139    *  WARNING: spatial information for SSD/DIGITS is obtained by pairing
 140
 141    *          digits on p and n side originating from the same track, when
 142
 143    *          possible. This (mis)use of DIGITS is tolerated for debugging
 144
 145    *          purposes only !!!!  The pairing in real life should be done
 146
 147    *          starting from info really available...
 148
 149    *
 150
 151    *  COMPILATION: this macro can be compiled.
 152
 153    *      1)       You need to set your include path with
 154
 155    * gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g");
 156
 157    *      3)       If you are using root instead of aliroot you need to
 158
 159    *               execute the macro loadlibs.C in advance
 160
 161    *      4)       To compile this macro from root (or aliroot):
 162
 163    *                 ---  .L AliITSGeoPlot.C++
 164
 165    *                 ---  AliITSGeoPlot();
 166
 167    *
 168
 169    *  M.Masera  14/05/2001 18:30
 170
 171    *  Last rev. 31/05/2004 14:00 (Clusters V2 added)  E.C.
 172
 173    ********************************************************************/
 174
 175
 176
 177   //Options
 178
 179   TString choice(opt);
 180
 181   Bool_t All = choice.Contains("All");
 182
 183   Bool_t verbose=choice.Contains("Verbose");
 184
 185   Bool_t userec=choice.Contains("Rec");
 186
 187   Bool_t useclustersv2=choice.Contains("ClustersV2");
 188
 189   Int_t retcode=1; //return code
 190
 191
 192
 193   if (gClassTable->GetID("AliRun") < 0) {
 194
 195     gInterpreter->ExecuteMacro("loadlibs.C");
 196
 197   }
 198
 199   else {
 200
 201     if(gAlice){
 202
 203       delete gAlice->GetRunLoader();
 204
 205       delete gAlice;
 206
 207       gAlice=0;
 208
 209     }
 210
 211   }
 212
 213
 214
 215   AliRunLoader* rl = AliRunLoader::Open(filename);
 216
 217   if (rl == 0x0){
 218
 219     cerr<<"AliITSGeoPlot.C : Can not open session RL=NULL"<< endl;
 220
 221     return -1;
 222
 223   }
 224
 225   Int_t retval = rl->LoadgAlice();
 226
 227   if (retval){
 228
 229     cerr<<"AliITSGeoPlot.C : LoadgAlice returned error"<<endl;
 230
 231     return -1;
 232
 233   }
 234
 235   gAlice=rl->GetAliRun();
 236
 237
 238
 239   retval = rl->LoadHeader();
 240
 241   if (retval){
 242
 243     cerr<<"AliITSGeoPlot.C : LoadHeader returned error"<<endl;
 244
 245     return -1;
 246
 247   }
 248
 249
 250
 251   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 252
 253
 254
 255   if(!ITSloader){
 256
 257     cerr<<"AliITSGeoPlot.C :  ITS loader not found"<<endl;
 258
 259     return -1;
 260
 261   }
 262
 263
 264
 265   ITSloader->LoadHits("read");
 266
 267   ITSloader->LoadDigits("read");
 268
 269   if(isfastpoints==1)ITSloader->SetRecPointsFileName("ITS.FastRecPoints.root");
 270
 271   ITSloader->LoadRecPoints("read");
 272
 273   rl->GetEvent(evesel);
 274
 275   Int_t nparticles = rl->GetHeader()->GetNtrack();
 276
 277   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 278
 279   ITS->SetTreeAddress();
 280
 281   if(verbose) {
 282
 283     cout<<" "<<endl<<" "<<endl;
 284
 285     cout<<"******* Event processing started   *******"<<endl;
 286
 287     cout<<"In the following, hits with 'StatusEntering' flag active"<<endl;
 288
 289     cout<<"will not be processed"<<endl;
 290
 291     cout << "Number of particles=  " << nparticles <<endl;
 292
 293   }
 294
 295
 296
 297   // HITS
 298
 299   TTree *TH = ITSloader->TreeH();
 300
 301   Stat_t ntracks = TH->GetEntries();
 302
 303   if(verbose)cout<<"Number of primary tracks= "<<ntracks<<endl;
 304
 305
 306
 307   // ITS
 308
 309   Int_t nmodules;
 310
 311   ITS->InitModules(-1,nmodules);
 312
 313   cout<<"Number of ITS modules= "<<nmodules<<endl;
 314
 315   cout<<"Filling modules... It takes a while, now. Please be patient"<<endl;
 316
 317   ITS->FillModules(0,0,nmodules," "," ");
 318
 319   cout<<"ITS modules .... DONE!"<<endl;
 320
 321
 322
 323   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec(ITSloader);
 324
 325   detTypeRec->SetDefaults();
 326
 327   // DIGITS
 328
 329   TTree *TD = ITSloader->TreeD();
 330
 331
 332
 333   //RECPOINTS (V2 clusters)
 334
 335   TTree *TR = ITSloader->TreeR();
 336
 337   TClonesArray *ITSrec  = detTypeRec->RecPoints();
 338
 339   TBranch *branch = 0;
 340
 341   if(userec && TR && ITSrec){
 342
 343     if(isfastpoints==1){
 344
 345       branch = ITSloader->TreeR()->GetBranch("ITSRecPointsF");
 346
 347       cout<<"using fast points\n";
 348
 349     }
 350
 351     else {
 352
 353       branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
 354
 355     }
 356
 357     if(branch)branch->SetAddress(&ITSrec);
 358
 359   }
 360
 361
 362
 363   if(userec && (!TR || !ITSrec || !branch)){
 364
 365     userec = kFALSE;
 366
 367     cout<<"\n ======================================================= \n";
 368
 369     cout<<"WARNING: there are no RECPOINTS on this file ! \n";
 370
 371     cout<<"======================================================= \n \n";
 372
 373   }
 374
 375   if(useclustersv2 && TR && ITSrec){
 376
 377     branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
 378
 379     if(branch)branch->SetAddress(&ITSrec);
 380
 381   }
 382
 383
 384
 385   if(useclustersv2 && (!TR || !ITSrec || !branch)){
 386
 387     useclustersv2 = kFALSE;
 388
 389     cout<<"\n ======================================================= \n";
 390
 391     cout<<"WARNING: there are no CLUSTERSV2 on this file ! \n";
 392
 393     cout<<"======================================================= \n \n";
 394
 395   }
 396
 397
 398
 399
 400
 401   //local variables
 402
 403   Int_t mod;   //module number
 404
 405   Int_t nbytes = 0;
 406
 407   Double_t pos[3];  // Global position of the current module
 408
 409   Float_t ragdet; // Radius of detector (x y plane)
 410
 411   Int_t first; // first module
 412
 413   Int_t last; // last module
 414
 415   Int_t nrecp; //number of RecPoints for a given module
 416
 417
 418
 419   //List of histograms
 420
 421   TObjArray histos(26,0);  // contains the pointers to the histograms
 422
 423   // Book histograms SPD
 424
 425   TH2F *rspd = new TH2F("rspd","Radii of digits - SPD",50,-10.,10.,50,-10.,10.);
 426
 427   TH2F *rhspd = new TH2F("rhspd","Radii of hits - SPD",50,-10.,10.,50,-10.,10.);
 428
 429   TH2F *rmspd = new TH2F("rmspd","Radii of SPD modules",50,-10.,10.,50,-10.,10.);
 430
 431   TH1F *zspd = new TH1F("zspd","Z of digits - SPD",100,-30.,30.);
 432
 433   TH1F *zhspd = new TH1F("zhspd","Z of hits - SPD",100,-30.,30.);
 434
 435   TH1F *zmspd = new TH1F("zmspd","Z of SPD modules",100,-30,30.);
 436
 437
 438
 439   Char_t title1[50]="";
 440
 441   Char_t title2[50]="";
 442
 443   if(userec){
 444
 445     sprintf(title1,"Radii of recpoints - %s","SPD");
 446
 447     sprintf(title2,"Z of recpoints - %s","SPD");
 448
 449   }
 450
 451   if(useclustersv2){
 452
 453     sprintf(title1,"Radii of clustersV2 - %s","SPD");
 454
 455     sprintf(title2,"Z of clustersV2 - %s","SPD");
 456
 457   }
 458
 459   TH2F *rrspd = new TH2F("rrspd",title1,50,-10.,10.,50,-10.,10.);
 460
 461   TH1F *zrspd = new TH1F("zrspd",title2,100,-30.,30.);
 462
 463   TH1F *enespd = new TH1F("enespd","Energy deposition SPD (KeV)",100,0.,1000.);
 464
 465   histos.AddLast(rspd);  // 0
 466
 467   histos.AddLast(rhspd); // 1
 468
 469   histos.AddLast(rmspd); // 2
 470
 471   histos.AddLast(zspd);  // 3
 472
 473   histos.AddLast(zhspd); // 4
 474
 475   histos.AddLast(zmspd); // 5
 476
 477   histos.AddLast(rrspd); // 6
 478
 479   histos.AddLast(zrspd); // 7
 480
 481   histos.AddLast(enespd); // 8
 482
 483   // Book histograms SDD
 484
 485   TH2F *rsdd = new TH2F("rsdd","Radii of digits - SDD",50,-40.,40.,50,-40.,40.);
 486
 487   TH2F *rhsdd = new TH2F("rhsdd","Radii of hits - SDD",50,-40.,40.,50,-40.,40.);
 488
 489   TH2F *rmsdd = new TH2F("rmsdd","Radii of SDD modules",50,-40.,40.,50,-40.,40.);
 490
 491   TH1F *zsdd = new TH1F("zsdd","Z of digits - SDD",100,-40.,40.);
 492
 493   TH1F *zhsdd = new TH1F("zhsdd","Z of hits - SDD",100,-40.,40.);
 494
 495   TH1F *zmsdd = new TH1F("zmsdd","Z of SDD modules",100,-40,40.);
 496
 497   Char_t title3[50];
 498
 499   Char_t title4[50];
 500
 501   if(userec){
 502
 503     sprintf(title3,"Radii of recpoints - %s","SDD");
 504
 505     sprintf(title4,"Z of recpoints - %s","SDD");
 506
 507   }
 508
 509   if(useclustersv2){
 510
 511     sprintf(title3,"Radii of clustersV2 - %s","SDD");
 512
 513     sprintf(title4,"Z of clustersV2 - %s","SDD");
 514
 515   }
 516
 517   TH2F *rrsdd = new TH2F("rrsdd",title3,50,-40.,40.,50,-40.,40.);
 518
 519   TH1F *zrsdd = new TH1F("zrsdd",title4,100,-40.,40.);
 520
 521   TH1F *enesdd = new TH1F("enesdd","Energy deposition SDD (KeV)",100,0.,1000.);
 522
 523   histos.AddLast(rsdd);  // 9
 524
 525   histos.AddLast(rhsdd); // 10
 526
 527   histos.AddLast(rmsdd); // 11
 528
 529   histos.AddLast(zsdd);  // 12
 530
 531   histos.AddLast(zhsdd); // 13
 532
 533   histos.AddLast(zmsdd); // 14
 534
 535   histos.AddLast(rrsdd); // 15
 536
 537   histos.AddLast(zrsdd); // 16
 538
 539   histos.AddLast(enesdd); // 17
 540
 541   // Book histogram SSD
 542
 543   TH2F *rssd = new TH2F("rssd","Radii of digits - SSD",50,-50.,50.,50,-50.,50.);
 544
 545   TH2F *rhssd = new TH2F("rhssd","Radii of hits - SSD",50,-50.,50.,50,-50.,50.);
 546
 547   TH2F *rmssd = new TH2F("rmssd","Radii of SSD modules",50,-50.,50.,50,-50.,50.);
 548
 549   TH1F *zssd = new TH1F("zssd","Z of digits - SSD",100,-70.,70.);
 550
 551   TH1F *zhssd = new TH1F("zhssd","Z of hits - SSD",100,-70.,70.);
 552
 553   TH1F *zmssd = new TH1F("zmssd","Z of SSD modules",100,-70,70.);
 554
 555   Char_t title5[50];
 556
 557   Char_t title6[50];
 558
 559   if(userec){
 560
 561     sprintf(title5,"Radii of recpoints - %s","SSD");
 562
 563     sprintf(title6,"Z of recpoints - %s","SSD");
 564
 565   }
 566
 567   if(useclustersv2){
 568
 569     sprintf(title5,"Radii of clustersV2 - %s","SSD");
 570
 571     sprintf(title6,"Z of clustersV2 - %s","SSD");
 572
 573   }
 574
 575
 576
 577   TH2F *rrssd = new TH2F("rrssd",title5,50,-50.,50.,50,-50.,50.);
 578
 579   TH1F *zrssd = new TH1F("zrssd",title6,100,-70.,70.);
 580
 581   TH1F *enessd = new TH1F("enessd","Energy deposition SSD (KeV)",100,0.,1000.);
 582
 583   histos.AddLast(rssd);  // 18
 584
 585   histos.AddLast(rhssd); // 19
 586
 587   histos.AddLast(rmssd); // 20
 588
 589   histos.AddLast(zssd);  // 21
 590
 591   histos.AddLast(zhssd); // 22
 592
 593   histos.AddLast(zmssd); // 23
 594
 595   histos.AddLast(rrssd); // 24
 596
 597   histos.AddLast(zrssd); // 25
 598
 599   histos.AddLast(enessd); // 26
 600
 601   //
 602
 603   // Loop on subdetectors
 604
 605   //
 606
 607   AliITSgeom *geom = ITS->GetITSgeom();
 608
 609   TString detna; // subdetector name
 610
 611   for(Int_t subd=0;subd<3;subd++){
 612
 613     if(All || (choice.Contains("SPD") && subd==0) || (choice.Contains("SDD") && subd==1) || (choice.Contains("SSD") && subd==2)){
 614
 615       // Prepare array for the digits
 616
 617       TClonesArray *ITSdigits  = ITS->DigitsAddress(subd);
 618
 619       Bool_t usedigits = kTRUE;
 620
 621       if(!ITSdigits){
 622
 623         usedigits = kFALSE;
 624
 625         cout<<"\n ======================================================= \n";
 626
 627         cout<<"WARNING: there are no DIGITS on this file ! \n";
 628
 629         cout<<"======================================================= \n \n";
 630
 631       }
 632
 633       // Get segmentation model
 634
 635       if(subd==0)detna="SPD";
 636
 637       if(subd==1)detna="SDD";
 638
 639       if(subd==2)detna="SSD";
 640
 641       AliITSsegmentation *seg=(AliITSsegmentation*)detTypeRec->GetSegmentationModel(subd);
 642
 643       // Loop on modules
 644
 645       first = geom->GetStartDet(subd);
 646
 647       last = geom->GetLastDet(subd);
 648
 649       if(verbose){
 650
 651         cout<<"     "<<endl<<"-------------------------------------"<<endl;
 652
 653         cout<<"Start processing subdetector "<<detna<<endl;
 654
 655         cout<<detna<<" first module "<<first<<endl;
 656
 657         cout<<detna<<" last module "<<last<<endl;
 658
 659         cout<<" "<<endl<<" "<<endl;
 660
 661       }
 662
 663       for (mod=first; mod<=last; mod++){
 664
 665         geom->GetTrans(mod,pos);  // position of the module in the MRS
 666
 667         ragdet=sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
 668
 669         // The following 2 histos are a check of the geometry
 670
 671         TH2F *bidi = (TH2F*)histos.At(2+subd*9);
 672
 673         TH1F *uni = (TH1F*)histos.At(5+subd*9);
 674
 675         bidi->Fill(pos[0],pos[1]);
 676
 677         uni->Fill(pos[2]);
 678
 679         if(verbose){
 680
 681           cout<<"=========================================================\n";
 682
 683           cout<<detna<<" module="<<mod<<endl;
 684
 685           cout<<"Mod. coordinates: "<<pos[0]<<", "<<pos[1]<<", ";
 686
 687           cout<<pos[2]<<" Radius "<<ragdet<<endl;
 688
 689         }
 690
 691
 692
 693         // Hits
 694
 695         GetHitsCoor(ITS,mod,histos,subd,verbose);
 696
 697
 698
 699         //RecPoints
 700
 701         if(userec){
 702
 703           detTypeRec->ResetRecPoints();
 704
 705           branch->GetEvent(mod);
 706
 707           TH2F *bidi=(TH2F*)histos.At(6+subd*9);
 708
 709           TH1F *uni=(TH1F*)histos.At(7+subd*9);
 710
 711           nrecp=GetRecCoor(geom,ITSrec,mod,bidi,uni,verbose);
 712
 713         }
 714
 715         if(useclustersv2){
 716
 717           detTypeRec->ResetRecPoints();
 718
 719           branch->GetEvent(mod);
 720
 721           TH2F *bidi=(TH2F*)histos.At(6+subd*9);
 722
 723           TH1F *uni=(TH1F*)histos.At(7+subd*9);
 724
 725           nrecp=GetClusCoor(geom,ITSrec,mod,bidi,uni,verbose);
 726
 727
 728
 729         }
 730
 731
 732
 733         // Digits
 734
 735         if(usedigits){
 736
 737           detTypeRec->ResetDigits();
 738
 739           nbytes += TD->GetEvent(mod);
 740
 741           GetDigits(seg,geom,ITSdigits,subd,mod,verbose,histos);
 742
 743         }
 744
 745
 746
 747       } // End of loop on the modules
 748
 749       TH1F *h1tmp;
 750
 751       TH2F *h2tmp;
 752
 753       //  Plot the histograms
 754
 755       TCanvas *current=0; // current Canvas (1--> SPD, 2---> SDD, 3---> SSD)
 756
 757       if(subd==0){
 758
 759         // Prepare canvas 1
 760
 761         TCanvas *c1 = new TCanvas("c1","SPD",10,10,600,900);
 762
 763         c1->Divide(2,3);
 764
 765         current=c1;
 766
 767       }
 768
 769       if(subd==1){
 770
 771         // Prepare canvas 2
 772
 773         TCanvas *c2 = new TCanvas("c2","SDD",40,40,600,900);
 774
 775         c2->Divide(2,3);
 776
 777         current=c2;
 778
 779       }
 780
 781       if(subd==2){
 782
 783         // Prepare canvas 3
 784
 785         TCanvas *c3 = new TCanvas("c3","SSD",70,70,600,900);
 786
 787         c3->Divide(2,3);
 788
 789         current=c3;
 790
 791       }
 792
 793       current->cd(1);
 794
 795       h2tmp = (TH2F*)histos.At(9*subd);
 796
 797       h2tmp->Draw();
 798
 799       current->cd(2);
 800
 801       h1tmp=(TH1F*)histos.At(3+subd*9);
 802
 803       h1tmp->Draw();
 804
 805       current->cd(3);
 806
 807       h2tmp=(TH2F*)histos.At(1+9*subd);
 808
 809       h2tmp->Draw();
 810
 811       current->cd(4);
 812
 813       h1tmp=(TH1F*)histos.At(4+subd*9);
 814
 815       h1tmp->Draw();
 816
 817
 818
 819       if(userec || useclustersv2){
 820
 821         current->cd(5);
 822
 823         h2tmp=(TH2F*)histos.At(6+9*subd);
 824
 825         h2tmp->Draw();
 826
 827         current->cd(6);
 828
 829         h1tmp=(TH1F*)histos.At(7+subd*9);
 830
 831         h1tmp->Draw();
 832
 833       }
 834
 835
 836
 837       else {
 838
 839         current->cd(5);
 840
 841         h2tmp=(TH2F*)histos.At(2+9*subd);
 842
 843         h2tmp->Draw();
 844
 845         current->cd(6);
 846
 847         h2tmp=(TH2F*)histos.At(5+9*subd);
 848
 849         h2tmp->Draw();
 850
 851       }
 852
 853     } // if(All.....
 854
 855   } // end of loop on subdetectors
 856
 857   // Save the histograms
 858
 859   TFile *fh = new TFile("AliITSGeoPlot.root","recreate");
 860
 861   // The list is written to file as a single entry
 862
 863   TList *lihi = new TList();
 864
 865   // copy the pointers to the histograms to a TList object.
 866
 867   // The histograms concerning recpoints are not copied if
 868
 869   // 'userec' is false.
 870
 871   for(Int_t i=0;i<histos.GetEntriesFast();i++){
 872
 873     if(choice.Contains("All") || (choice.Contains("SPD") && i<8) || (choice.Contains("SDD") && i>7 && i<16) || (choice.Contains("SSD") && i>15)){
 874
 875       if(!(!userec && ((i+2)%9==0 || (i+1)%9==0)))lihi->Add(histos.At(i));
 876
 877     }
 878
 879   }
 880
 881   lihi->Write("Histograms ITS hits+digits+recpoints",TObject::kSingleKey);
 882
 883   fh->Close();
 884
 885
 886
 887   return retcode;
 888
 889 }
 890
 891
 892
 893
 894
 895 void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb){
 896
 897   TH2F *h2=(TH2F*)histos.At(1+subd*9);
 898
 899   TH1F *h1=(TH1F*)histos.At(4+subd*9);
 900
 901   TH1F *ener=(TH1F*)histos.At(8+subd*9);
 902
 903   AliITS *ITS= (AliITS*)its;
 904
 905   AliITSmodule *modu = ITS->GetModule(mod);
 906
 907   TObjArray *fHits = modu->GetHits();
 908
 909   Int_t nhits = fHits->GetEntriesFast();
 910
 911   if(nhits>0){
 912
 913     if(verb){
 914
 915       cout<<"-------------------------------------------------------"<<endl;
 916
 917       cout<<"Number of HITS for module "<<mod<<": "<<nhits<<endl;
 918
 919     }
 920
 921     for (Int_t hit=0;hit<nhits;hit++) {
 922
 923       AliITShit *iHit = (AliITShit*) fHits->At(hit);
 924
 925       if(!iHit->StatusEntering()){
 926
 927         Float_t x=iHit->GetXG();
 928
 929         Float_t y=iHit->GetYG();
 930
 931         Float_t z=iHit->GetZG();
 932
 933         Float_t edep=iHit->GetIonization()*1000000;
 934
 935         h2->Fill(x,y);
 936
 937         h1->Fill(z);
 938
 939         ener->Fill(edep);
 940
 941         if(verb){
 942
 943           cout<<"hit # "<<hit<<" Global coordinates "<<x<<" "<<y<<" "<<z<<endl;
 944
 945           cout<<"track # "<<iHit->GetTrack()<<" energy deposition (KeV)= ";
 946
 947           cout<<edep<<endl;
 948
 949         }
 950
 951       }
 952
 953     }
 954
 955   }
 956
 957 }
 958
 959
 960
 961
 962
 963 Int_t GetClusCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
 964
 965
 966
 967   AliITSgeom *geom = (AliITSgeom*)ge;
 968
 969   Int_t nrecp = ITSrec->GetEntries();
 970
 971   if(nrecp>0){
 972
 973     Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
 974
 975     Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
 976
 977     if(verb){
 978
 979       cout<<"-------------------------------------------------------"<<endl;
 980
 981       cout<<"Number of CLUSTERS for module "<<mod<<": "<<nrecp<<endl;
 982
 983     }
 984
 985     for(Int_t irec=0;irec<nrecp;irec++) {
 986
 987       AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
 988
 989       Double_t rot[9];
 990
 991       geom->GetRotMatrix(mod,rot);
 992
 993       Int_t lay,lad,det;
 994
 995       geom->GetModuleId(mod,lay,lad,det);
 996
 997       Float_t tx,ty,tz;
 998
 999       geom->GetTrans(lay,lad,det,tx,ty,tz);
1000
1001
1002
1003       Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1004
1005       Double_t phi1=TMath::Pi()/2+alpha;
1006
1007       if(lay==1) phi1+=TMath::Pi();
1008
1009
1010
1011       Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1012
1013       Float_t  r=tx*cp+ty*sp;
1014
1015       gc[0]= r*cp - recp->GetY()*sp;
1016
1017       gc[1]= r*sp + recp->GetY()*cp;
1018
1019       gc[2]=recp->GetZ();
1020
1021
1022
1023       if(verb){
1024
1025         Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
1026
1027         cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" ";
1028
1029         cout<<r<<endl;
1030
1031         cout<<"Associated track "<<recp->GetLabel(0)<<endl;
1032
1033       }
1034
1035       h2->Fill(gc[0],gc[1]);
1036
1037       h1->Fill(gc[2]);
1038
1039
1040
1041     }
1042
1043   }
1044
1045   return nrecp;
1046
1047 }
1048
1049 Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
1050
1051   AliITSgeom *geom = (AliITSgeom*)ge;
1052
1053   Int_t nrecp = ITSrec->GetEntries();
1054
1055   if(nrecp>0){
1056
1057     Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
1058
1059     Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
1060
1061     if(verb){
1062
1063       cout<<"-------------------------------------------------------"<<endl;
1064
1065       cout<<"Number of REC POINTS for module "<<mod<<": "<<nrecp<<endl;
1066
1067     }
1068
1069     for(Int_t irec=0;irec<nrecp;irec++) {
1070
1071       AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
1072
1073       lc[0]=recp->GetDetLocalX();
1074
1075       lc[2]=recp->GetDetLocalZ();
1076
1077       geom->LtoG(mod,lc,gc);
1078
1079       if(verb){
1080
1081         cout<<"recp # "<<irec<<" local coordinates. lx= "<<lc[0]<<" lz= ";
1082
1083         cout<<lc[2]<<endl;
1084
1085         Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
1086
1087         cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" ";
1088
1089         cout<<r<<endl;
1090
1091         cout<<"Associated track "<<recp->GetLabel(0)<<endl;
1092
1093       }
1094
1095       h2->Fill(gc[0],gc[1]);
1096
1097       h1->Fill(gc[2]);
1098
1099     }
1100
1101   }
1102
1103   return nrecp;
1104
1105 }
1106
1107
1108
1109 void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos){
1110
1111   AliITSsegmentation *seg = (AliITSsegmentation*)tmps;
1112
1113   AliITSgeom *geom = (AliITSgeom*)ge;
1114
1115   Int_t layer;
1116
1117   Int_t ladder;
1118
1119   Int_t detec;
1120
1121   if(subd==2){
1122
1123     geom->GetModuleId(mod,layer,ladder,detec);
1124
1125     seg->SetLayer(layer);
1126
1127   }
1128
1129   Float_t lcoor[3]; for(Int_t j=0; j<3; j++) lcoor[j]=0.;  //local coord dig.
1130
1131   Float_t gcoor[3]; for(Int_t j=0; j<3; j++) gcoor[j]=0.; // global coo. dig.
1132
1133   Float_t ragdig; // Radius digit
1134
1135   TArrayI ssdone(5000);  // used to match p and n side digits of strips
1136
1137   TArrayI pair(5000);    // as above
1138
1139   Int_t ndigits = ITSdigits->GetEntries();
1140
1141   AliITSdigit *digs;
1142
1143   if(ndigits){
1144
1145     if(verbose){
1146
1147       cout<<"-------------------------------------------------------"<<endl;
1148
1149       cout<<"Number of DIGITS for module "<<mod<<": "<<ndigits<<endl;
1150
1151     }
1152
1153     // Get the coordinates of the module
1154
1155     if(subd==2){
1156
1157       for (Int_t digit=0;digit<ndigits;digit++){
1158
1159         ssdone[digit]=0;
1160
1161         pair[digit]=0;
1162
1163       }
1164
1165     }
1166
1167     for (Int_t digit=0;digit<ndigits;digit++) {
1168
1169       digs = (AliITSdigit*)ITSdigits->UncheckedAt(digit);
1170
1171       Int_t iz=digs->GetCoord1();  // cell number z
1172
1173       Int_t ix=digs->GetCoord2();  // cell number x
1174
1175       // Get local coordinates of the element
1176
1177       if(subd<2){
1178
1179         seg->DetToLocal(ix,iz,lcoor[0],lcoor[2]);
1180
1181       }
1182
1183       else{
1184
1185         // SSD: if iz==0 ---> N side; if iz==1 P side
1186
1187         if(ssdone[digit]==0){
1188
1189           ssdone[digit]=1;
1190
1191           pair[digit]=-1;
1192
1193           Bool_t pside=(iz==1);
1194
1195           Bool_t impaired=kTRUE;
1196
1197           Int_t pstrip=0;
1198
1199           Int_t nstrip=0;
1200
1201           if(pside)pstrip=ix;
1202
1203           if(!pside)nstrip=ix;
1204
1205           for(Int_t digi2=0;digi2<ndigits;digi2++){
1206
1207             if(ssdone[digi2]==0 && impaired){
1208
1209               AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(digi2);
1210
1211               if(dig2->GetCoord1() != iz && dig2->GetTrack(0)==digs->GetTrack(0) && dig2->GetTrack(0)>0){
1212
1213                 ssdone[digi2]=2;
1214
1215                 pair[digit]=digi2;
1216
1217                 if(pside)nstrip=dig2->GetCoord2();
1218
1219                 if(!pside)pstrip=dig2->GetCoord2();
1220
1221                 impaired=kFALSE;
1222
1223               }
1224
1225             }
1226
1227           }
1228
1229           if(!impaired)seg->GetPadCxz(pstrip,nstrip,lcoor[0],lcoor[2]);
1230
1231         }
1232
1233       }
1234
1235       if(subd<2 || (subd==2 && ssdone[digit]==1)){
1236
1237         Int_t coor1=digs->GetCoord1();
1238
1239         Int_t coor2=digs->GetCoord2();
1240
1241         Int_t tra0=digs->GetTrack(0);
1242
1243         if(verbose){
1244
1245           cout<<"digit # "<<digit<<" fCoord1= "<<coor1<<" fCoord2= ";
1246
1247           cout<<coor2<<" track "<<tra0<<" "<<endl;
1248
1249           if(subd<2)cout<<"local coordinates -- x="<<lcoor[0]<<", z=";
1250
1251           cout<<lcoor[2]<<endl;
1252
1253           if(subd==2){
1254
1255             if(pair[digit]==-1){
1256
1257               cout<<"(digit  not paired)"<<endl;
1258
1259             }
1260
1261             else {
1262
1263               cout<<"local coordinates -- x="<<lcoor[0]<<", z="<<lcoor[2];
1264
1265               cout<<endl;
1266
1267               Int_t dtmp=pair[digit];
1268
1269               AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(dtmp);
1270
1271               Int_t coor1b=dig2->GetCoord1();
1272
1273               Int_t coor2b=dig2->GetCoord2();
1274
1275               Int_t tra0b=dig2->GetTrack(0);
1276
1277               cout<<"(digit paired with digit #"<<dtmp<<endl;
1278
1279               cout<<"with fCoord1= "<<coor1b<<" fCoord2= "<<coor2b;
1280
1281               cout<<" track "<<tra0b<<")"<<endl;
1282
1283             }
1284
1285           }
1286
1287         }
1288
1289         if(subd<2 || (subd==2 && pair[digit]!=-1)){
1290
1291           // Global coordinates of the element
1292
1293           //SDD and SPD use cm, SSD microns (GetPadCxz)
1294
1295           if(subd==2)for(Int_t j=0;j<3;j++)lcoor[j]=lcoor[j]/10000.;
1296
1297           lcoor[1]=0.;
1298
1299           geom->LtoG(mod,lcoor,gcoor);  // global coord. in cm
1300
1301           ragdig=sqrt(gcoor[0]*gcoor[0]+gcoor[1]*gcoor[1]);
1302
1303           if(verbose){
1304
1305             cout<<"global coordinates "<<gcoor[0]<<" "<<gcoor[1];
1306
1307             cout<<" "<<gcoor[2]<<" Radius "<<ragdig<<endl;
1308
1309           }
1310
1311           //Fill histograms
1312
1313           TH2F *bidi = (TH2F*)histos.At(subd*9);
1314
1315           TH1F *uni = (TH1F*)histos.At(3+subd*9);
1316
1317           bidi->Fill(gcoor[0],gcoor[1]);
1318
1319           uni->Fill(gcoor[2]);
1320
1321         }
1322
1323       }
1324
1325     } // loop on digits for this module
1326
1327   } // if(ndigits>0....
1328
1329 }
1330