ITS/AliITS.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.49  2001/05/11 07:37:49  hristov
  19 Legacy lines commented
  20
  21 Revision 1.48  2001/05/10 18:14:25  barbera
  22 A typo corrected
  23
  24 Revision 1.47  2001/05/10 17:55:59  barbera
  25 Modified to create rec points also for PPR geometries
  26
  27 Revision 1.46  2001/05/10 00:05:28  nilsen
  28 Allowed for HitsToDigits function to work with versions 5, 7, 8, and 9. This
  29 should probably be cleaned up to only check to make sure that fITSgeom has
  30 been properly defined.
  31
  32 Revision 1.45  2001/05/01 22:35:48  nilsen
  33 Remove/commented a number of cout<< statements. and made change needed by
  34 SSD code.
  35
  36 Revision 1.44  2001/04/26 22:44:01  nilsen
  37 Removed dependence on layer 5/6 in AliITS::HitsToDigits. This will be
  38 done properly in AliITSv???.cxx via SetDefaults.
  39
  40 Revision 1.43  2001/04/26 13:22:52  barbera
  41 TMatrix and TVector elimininated to speed up the code
  42
  43 Revision 1.42  2001/04/25 21:55:12  barbera
  44 Updated version to be compatible with actual verion of STEER and TPC
  45
  46 Revision 1.41  2001/04/21 15:16:51  barbera
  47 Updated with the new SSD reconstruction code
  48
  49 Revision 1.40  2001/03/17 15:07:06  mariana
  50 Update SDD response parameters
  51
  52 Revision 1.39  2001/03/12 17:45:32  hristov
  53 Changes needed on Sun with CC 5.0
  54
  55 Revision 1.38  2001/03/07 14:04:51  barbera
  56 Some vector dimensions increased to cope with full events
  57
  58 Revision 1.37  2001/03/07 12:36:35  barbera
  59 A change added in the tracking part to manage delta rays
  60
  61 Revision 1.36  2001/03/02 19:44:11  barbera
  62  modified to taking into account new version tracking v1
  63
  64 Revision 1.35  2001/02/28 18:16:46  mariana
  65 Make the code compatible with the new AliRun
  66
  67 Revision 1.34  2001/02/11 15:51:39  mariana
  68 Set protection in MakeBranch
  69
  70 Revision 1.33  2001/02/10 22:26:39  mariana
  71 Move the initialization of the containers for raw clusters in MakeTreeC()
  72
  73 Revision 1.32  2001/02/08 23:55:31  nilsen
  74 Removed fMajor/MinorVersion variables in favor of variables in derived classes.
  75 Set arrays char *det[3] = {"SPD","SDD","SSD"} as const.
  76
  77 Revision 1.31  2001/02/02 23:57:28  nilsen
  78 Added include file that are no londer included in AliITSgeom.h
  79
  80 Revision 1.30  2001/01/30 09:23:13  hristov
  81 Streamers removed (R.Brun)
  82
  83 Revision 1.29  2001/01/26 20:01:09  hristov
  84 Major upgrade of AliRoot code
  85
  86 Revision 1.28  2000/12/18 14:02:00  barbera
  87 new version of the ITS tracking to take into account the new TPC track parametrization
  88
  89 Revision 1.27  2000/12/08 13:49:27  barbera
  90 Hidden declaration in a for loop removed to be compliant with HP-UX compiler
  91
  92 Revision 1.26  2000/11/27 13:12:13  barbera
  93 New version containing the files for tracking
  94
  95 Revision 1.25  2000/11/12 22:38:05  barbera
  96 Added header file for the SPD Bari model
  97
  98 Revision 1.24  2000/10/09 22:18:12  barbera
  99 Bug fixes from MAriana to le AliITStest.C run correctly
 100
 101 Revision 1.23  2000/10/05 20:47:42  nilsen
 102 fixed dependencies of include files. Tryed but failed to get a root automaticly
 103 generates streamer function to work. Modified SetDefaults.
 104
 105 Revision 1.9.2.15  2000/10/04 16:56:40  nilsen
 106 Needed to include stdlib.h
 107
 108 =======
 109 Revision 1.22  2000/10/04 19:45:52  barbera
 110 Corrected by F. Carminati for v3.04
 111
 112 Revision 1.21  2000/10/02 21:28:08  fca
 113 Removal of useless dependecies via forward declarations
 114
 115 Revision 1.20  2000/10/02 16:31:39  barbera
 116 General code clean-up
 117
 118 Revision 1.9.2.14  2000/10/02 15:43:51  barbera
 119 General code clean-up (e.g., printf -> cout)
 120
 121 Revision 1.19  2000/09/22 12:13:25  nilsen
 122 Patches and updates for fixes to this and other routines.
 123
 124 Revision 1.18  2000/07/12 05:32:20  fca
 125 Correcting several syntax problem with static members
 126
 127 Revision 1.17  2000/07/10 16:07:18  fca
 128 Release version of ITS code
 129
 130 Revision 1.9.2.3  2000/02/02 13:42:09  barbera
 131 fixed AliITS.cxx for new AliRun structure. Added ITS hits list to list of hits which will have their track numbers updated
 132
 133 Revision 1.9.2.2  2000/01/23 03:03:13  nilsen
 134 //fixed FillModule. Removed fi(fabs(xl)<dx....
 135
 136 Revision 1.9.2.1  2000/01/12 19:03:32  nilsen
 137 This is the version of the files after the merging done in December 1999.
 138 See the ReadMe110100.txt file for details
 139
 140 Revision 1.9  1999/11/14 14:33:25  fca
 141 Correct problems with distructors and pointers, thanks to I.Hrivnacova
 142
 143 Revision 1.8  1999/09/29 09:24:19  fca
 144 Introduction of the Copyright and cvs Log
 145
 146 */
 147
 148 ///////////////////////////////////////////////////////////////////////////////
 149 //
 150 //      An overview of the basic philosophy of the ITS code development
 151 // and analysis is show in the figure below.
 152 //Begin_Html
 153 /*
 154 <img src="picts/ITS/ITS_Analysis_schema.gif">
 155 </pre>
 156 <br clear=left>
 157 <font size=+2 color=red>
 158 <p>Roberto Barbera is in charge of the ITS Offline code (1999).
 159 <a href="mailto:roberto.barbera@ct.infn.it">Roberto Barbera</a>.
 160 </font>
 161 <pre>
 162 */
 163 //End_Html
 164 //
 165 //  AliITS. Inner Traking System base class.
 166 //  This class contains the base procedures for the Inner Tracking System
 167 //
 168 //Begin_Html
 169 /*
 170 <img src="picts/ITS/AliITS_Class_Diagram.gif">
 171 </pre>
 172 <br clear=left>
 173 <font size=+2 color=red>
 174 <p>This show the class diagram of the different elements that are part of
 175 the AliITS class.
 176 </font>
 177 <pre>
 178 */
 179 //End_Html
 180 //
 181 // Version: 0
 182 // Written by Rene Brun, Federico Carminati, and Roberto Barbera
 183 //
 184 // Version: 1
 185 // Modified and documented by Bjorn S. Nilsen
 186 // July 11 1999
 187 //
 188 // Version: 2
 189 // Modified and documented by A. Bologna
 190 // October 18 1999
 191 //
 192 // AliITS is the general base class for the ITS. Also see AliDetector for
 193 // futher information.
 194 //
 195 ///////////////////////////////////////////////////////////////////////////////
 196 #include <iostream.h>
 197 #include <iomanip.h>
 198 #include <fstream.h>
 199 #include <stdlib.h>
 200 #include <TMath.h>
 201 #include <TRandom.h>
 202 #include <TBranch.h>
 203 #include <TVector.h>
 204 #include <TClonesArray.h>
 205 #include <TROOT.h>
 206 #include <TObjectTable.h>
 207 #include <TFile.h>
 208 #include <TTree.h>
 209 #include <TString.h>
 210 #include <TParticle.h>
 211
 212
 213 #include "AliRun.h"
 214 #include "AliITS.h"
 215 #include "AliITSMap.h"
 216 #include "AliITSDetType.h"
 217 #include "AliITSClusterFinder.h"
 218 //#include "AliITSsimulation.h"
 219 #include "AliITSsimulationSPD.h"
 220 #include "AliITSsimulationSDD.h"
 221 #include "AliITSsimulationSSD.h"
 222 #include "AliITSresponse.h"
 223 #include "AliITSsegmentationSPD.h"
 224 #include "AliITSresponseSPD.h"
 225 #include "AliITSresponseSPDbari.h"
 226 #include "AliITSsegmentationSDD.h"
 227 #include "AliITSresponseSDD.h"
 228 #include "AliITSsegmentationSSD.h"
 229 #include "AliITSresponseSSD.h"
 230 #include "AliITShit.h"
 231 #include "AliITSgeom.h"
 232 #include "AliITSdigit.h"
 233 #include "AliITSmodule.h"
 234 #include "AliITSRecPoint.h"
 235 #include "AliITSRawCluster.h"
 236 #include "AliMC.h"
 237 #include "stdlib.h"
 238 #include "AliKalmanTrack.h"
 239 #include "AliMagF.h"
 240
 241 #include "AliITStrack.h"
 242 #include "AliITSiotrack.h"
 243 #include "AliITStracking.h"
 244 #include "AliITSRad.h"
 245 #include "../TPC/AliTPC.h"
 246 #include "../TPC/AliTPCParam.h"
 247
 248
 249 ClassImp(AliITS)
 250
 251 //_____________________________________________________________________________
 252 AliITS::AliITS() : AliDetector() {
 253   //
 254   // Default initialiser for ITS
 255   //     The default constructor of the AliITS class. In addition to
 256   // creating the AliITS class it zeros the variables fIshunt (a member
 257   // of AliDetector class), fEuclidOut, and fIdN, and zeros the pointers
 258   // fITSpoints, fIdSens, and fIdName. The AliDetector default constructor
 259   // is also called.
 260   //
 261
 262
 263   fIshunt     = 0;
 264   fEuclidOut  = 0;
 265
 266   fNDetTypes = kNTYPES;
 267   fIdN        = 0;
 268   fIdName     = 0;
 269   fIdSens     = 0;
 270   fITSmodules = 0;
 271   //
 272   fDetTypes   = 0;
 273   //
 274   fDtype  = 0;
 275   fNdtype = 0;
 276   fCtype  = 0;
 277   fNctype = 0;
 278   fRecPoints = 0;
 279   fNRecPoints = 0;
 280   fTreeC = 0;
 281   //
 282   fITSgeom=0;
 283 }
 284
 285 //_____________________________________________________________________________
 286 AliITS::AliITS(const char *name, const char *title):AliDetector(name,title){
 287   //
 288   // Default initialiser for ITS
 289   //     The constructor of the AliITS class. In addition to creating the
 290   // AliITS class, it allocates memory for the TClonesArrays fHits and
 291   // fDigits, and for the TObjArray fITSpoints. It also zeros the variables
 292   // fIshunt (a member of AliDetector class), fEuclidOut, and fIdN, and zeros
 293   // the pointers fIdSens and fIdName. To help in displaying hits via the ROOT
 294   // macro display.C AliITS also sets the marker color to red. The variables
 295   // passes with this constructor, const char *name and *title, are used by
 296   // the constructor of AliDetector class. See AliDetector class for a
 297   // description of these parameters and its constructor functions.
 298   //
 299
 300
 301   fHits       = new TClonesArray("AliITShit", 1560);
 302   gAlice->AddHitList(fHits);
 303
 304   fNDetTypes = kNTYPES;
 305
 306   fNdtype = new Int_t[kNTYPES];
 307   fDtype = new TObjArray(kNTYPES);
 308
 309   fNctype = new Int_t[kNTYPES];
 310   fCtype = new TObjArray(kNTYPES);
 311
 312
 313   fRecPoints = 0;
 314   fNRecPoints = 0;
 315
 316   fTreeC = 0;
 317
 318   fITSmodules = 0;
 319
 320   fIshunt     = 0;
 321   fEuclidOut  = 0;
 322   fIdN        = 0;
 323   fIdName     = 0;
 324   fIdSens     = 0;
 325
 326   fDetTypes = new TObjArray(kNTYPES);
 327
 328   Int_t i;
 329   for(i=0;i<kNTYPES;i++) {
 330     (*fDetTypes)[i]=new AliITSDetType();
 331     fNdtype[i]=0;
 332     fNctype[i]=0;
 333    }
 334   //
 335
 336   SetMarkerColor(kRed);
 337
 338   fITSgeom=0;
 339 }
 340 //___________________________________________________________________________
 341 AliITS::AliITS(AliITS &source){
 342   // copy constructor
 343   if(this==&source) return;
 344   Error("AliITS::Copy constructor",
 345         "You are not allowed to make a copy of the AliITS");
 346   exit(1);
 347 }
 348 //____________________________________________________________________________
 349 AliITS& AliITS::operator=(AliITS &source){
 350   // assignment operator
 351   if(this==&source) return *this;
 352   Error("AliITS::operator=",
 353         "You are not allowed to make a copy of the AliITS");
 354   exit(1);
 355   return *this; //fake return
 356 }
 357 //____________________________________________________________________________
 358 void AliITS::ClearModules(){
 359   //clear the modules TObjArray
 360
 361   if(fITSmodules) fITSmodules->Delete();
 362
 363 }
 364 //_____________________________________________________________________________
 365 AliITS::~AliITS(){
 366   //
 367   // Default distructor for ITS
 368   //     The default destructor of the AliITS class. In addition to deleting
 369   // the AliITS class it deletes the memory pointed to by the fHits, fDigits,
 370   // fIdSens, fIdName, and fITSpoints.
 371   //
 372
 373
 374   delete fHits;
 375   delete fDigits;
 376   delete fRecPoints;
 377 //  delete fIdName;        // TObjArray of TObjStrings
 378   if(fIdName!=0) delete[] fIdName;  // Array of TStrings
 379   if(fIdSens!=0) delete[] fIdSens;
 380   if(fITSmodules!=0) {
 381       this->ClearModules();
 382       delete fITSmodules;
 383   }// end if fITSmodules!=0
 384
 385   //
 386   if(fDtype) {
 387     fDtype->Delete();
 388     delete fDtype;
 389   }
 390   delete [] fNdtype;
 391   if (fCtype) {
 392     fCtype->Delete();
 393     delete fCtype;
 394   }
 395   delete [] fNctype;
 396   //
 397
 398   if (fDetTypes) {
 399     fDetTypes->Delete();
 400     delete fDetTypes;
 401   }
 402
 403   if (fTreeC) delete fTreeC;
 404
 405   if (fITSgeom) delete fITSgeom;
 406
 407 }
 408
 409 //___________________________________________
 410 AliITSDetType* AliITS::DetType(Int_t id)
 411 {
 412   //return pointer to id detector type
 413     return ((AliITSDetType*) (*fDetTypes)[id]);
 414
 415 }
 416 //___________________________________________
 417 void AliITS::SetClasses(Int_t id, const char *digit, const char *cluster)
 418 {
 419   //set the digit and cluster classes to be used for the id detector type
 420     ((AliITSDetType*) (*fDetTypes)[id])->ClassNames(digit,cluster);
 421
 422 }
 423 //___________________________________________
 424 void AliITS::SetResponseModel(Int_t id, AliITSresponse *response)
 425 {
 426   //set the response model for the id detector type
 427
 428     ((AliITSDetType*) (*fDetTypes)[id])->ResponseModel(response);
 429
 430 }
 431
 432 //___________________________________________
 433 void AliITS::SetSegmentationModel(Int_t id, AliITSsegmentation *seg)
 434 {
 435   //set the segmentation model for the id detector type
 436
 437     ((AliITSDetType*) (*fDetTypes)[id])->SegmentationModel(seg);
 438
 439 }
 440
 441 //___________________________________________
 442 void AliITS::SetSimulationModel(Int_t id, AliITSsimulation *sim)
 443 {
 444   //set the simulation model for the id detector type
 445
 446    ((AliITSDetType*) (*fDetTypes)[id])->SimulationModel(sim);
 447
 448 }
 449 //___________________________________________
 450 void AliITS::SetReconstructionModel(Int_t id, AliITSClusterFinder *reconst)
 451 {
 452   //set the cluster finder model for the id detector type
 453
 454    ((AliITSDetType*) (*fDetTypes)[id])->ReconstructionModel(reconst);
 455
 456 }
 457
 458 //_____________________________________________________________________________
 459 void AliITS::AddHit(Int_t track, Int_t *vol, Float_t *hits){
 460   //
 461   // Add an ITS hit
 462   //     The function to add information to the AliITShit class. See the
 463   // AliITShit class for a full description. This function allocates the
 464   // necessary new space for the hit information and passes the variable
 465   // track, and the pointers *vol and *hits to the AliITShit constructor
 466   // function.
 467   //
 468   TClonesArray &lhits = *fHits;
 469   new(lhits[fNhits++]) AliITShit(fIshunt,track,vol,hits);
 470 }
 471 //_____________________________________________________________________________
 472 void AliITS::AddRealDigit(Int_t id, Int_t *digits)
 473 {
 474   // add a real digit - as coming from data
 475
 476   TClonesArray &ldigits = *((TClonesArray*)(*fDtype)[id]);
 477   new(ldigits[fNdtype[id]++]) AliITSdigit(digits);
 478
 479 }
 480 //_____________________________________________________________________________
 481 void AliITS::AddSimDigit(Int_t id, AliITSdigit *d)
 482 {
 483
 484   // add a simulated digit
 485
 486   TClonesArray &ldigits = *((TClonesArray*)(*fDtype)[id]);
 487
 488   switch(id)
 489   {
 490   case 0:
 491      new(ldigits[fNdtype[id]++]) AliITSdigitSPD(*((AliITSdigitSPD*)d));
 492      break;
 493   case 1:
 494      new(ldigits[fNdtype[id]++]) AliITSdigitSDD(*((AliITSdigitSDD*)d));
 495      break;
 496   case 2:
 497      new(ldigits[fNdtype[id]++]) AliITSdigitSSD(*((AliITSdigitSSD*)d));
 498      break;
 499   }
 500
 501 }
 502
 503 //_____________________________________________________________________________
 504 void AliITS::AddSimDigit(Int_t id,Float_t phys,Int_t *digits,Int_t *tracks,Int_t *hits,Float_t *charges){
 505
 506   // add a simulated digit to the list
 507
 508   TClonesArray &ldigits = *((TClonesArray*)(*fDtype)[id]);
 509   switch(id)
 510   {
 511   case 0:
 512      new(ldigits[fNdtype[id]++]) AliITSdigitSPD(digits,tracks,hits);
 513      break;
 514   case 1:
 515      new(ldigits[fNdtype[id]++]) AliITSdigitSDD(phys,digits,tracks,hits,charges);
 516      break;
 517   case 2:
 518      new(ldigits[fNdtype[id]++]) AliITSdigitSSD(digits,tracks,hits);
 519      break;
 520   }
 521
 522 }
 523
 524 //_____________________________________________________________________________
 525 void AliITS::AddCluster(Int_t id, AliITSRawCluster *c)
 526 {
 527
 528   // add a cluster to the list
 529
 530   TClonesArray &lcl = *((TClonesArray*)(*fCtype)[id]);
 531
 532   switch(id)
 533   {
 534   case 0:
 535      new(lcl[fNctype[id]++]) AliITSRawClusterSPD(*((AliITSRawClusterSPD*)c));
 536      break;
 537   case 1:
 538      new(lcl[fNctype[id]++]) AliITSRawClusterSDD(*((AliITSRawClusterSDD*)c));
 539      break;
 540   case 2:
 541      new(lcl[fNctype[id]++]) AliITSRawClusterSSD(*((AliITSRawClusterSSD*)c));
 542      break;
 543   }
 544
 545 }
 546
 547
 548 //_____________________________________________________________________________
 549 void AliITS::AddRecPoint(const AliITSRecPoint &r)
 550 {
 551   //
 552   // Add a reconstructed space point to the list
 553   //
 554   TClonesArray &lrecp = *fRecPoints;
 555   new(lrecp[fNRecPoints++]) AliITSRecPoint(r);
 556 }
 557
 558
 559 //____________________________________________
 560 void AliITS::ResetDigits()
 561 {
 562     //
 563     // Reset number of digits and the digits array for the ITS detector
 564     //
 565
 566     if (!fDtype) return;
 567
 568     Int_t i;
 569     for (i=0;i<kNTYPES;i++ ) {
 570         if ((*fDtype)[i])    ((TClonesArray*)(*fDtype)[i])->Clear();
 571         if (fNdtype)  fNdtype[i]=0;
 572     }
 573 }
 574
 575 //____________________________________________
 576 void AliITS::ResetDigits(Int_t i)
 577 {
 578     //
 579     // Reset number of digits and the digits array for this branch
 580     //
 581   if ((*fDtype)[i])    ((TClonesArray*)(*fDtype)[i])->Clear();
 582   if (fNdtype)  fNdtype[i]=0;
 583 }
 584
 585
 586 //____________________________________________
 587 void AliITS::ResetClusters()
 588 {
 589     //
 590     // Reset number of clusters and the clusters array for ITS
 591     //
 592
 593     Int_t i;
 594     for (i=0;i<kNTYPES;i++ ) {
 595         if ((*fCtype)[i])    ((TClonesArray*)(*fCtype)[i])->Clear();
 596         if (fNctype)  fNctype[i]=0;
 597     }
 598
 599 }
 600
 601 //____________________________________________
 602 void AliITS::ResetClusters(Int_t i)
 603 {
 604     //
 605     // Reset number of clusters and the clusters array for this branch
 606     //
 607         if ((*fCtype)[i])    ((TClonesArray*)(*fCtype)[i])->Clear();
 608         if (fNctype)  fNctype[i]=0;
 609
 610 }
 611
 612
 613 //____________________________________________
 614 void AliITS::ResetRecPoints()
 615 {
 616     //
 617     // Reset number of rec points and the rec points array
 618     //
 619     if (fRecPoints) fRecPoints->Clear();
 620     fNRecPoints = 0;
 621
 622 }
 623
 624 //_____________________________________________________________________________
 625 Int_t AliITS::DistancetoPrimitive(Int_t , Int_t ){
 626   //
 627   // Distance from mouse to ITS on the screen. Dummy routine
 628   //     A dummy routine used by the ROOT macro display.C to allow for the
 629   // use of the mouse (pointing device) in the macro. In general this should
 630   // never be called. If it is it returns the number 9999 for any value of
 631   // x and y.
 632   //
 633   return 9999;
 634 }
 635
 636 //_____________________________________________________________________________
 637 void AliITS::Init(){
 638   //
 639   // Initialise ITS after it has been built
 640   //     This routine initializes the AliITS class. It is intended to be called
 641   // from the Init function in AliITSv?. Besides displaying a banner
 642   // indicating that it has been called it initializes the array fIdSens
 643   // and sets the default segmentation, response, digit and raw cluster classes
 644   // Therefore it should be called after a call to CreateGeometry.
 645   //
 646   Int_t i;
 647
 648 //
 649   SetDefaults();
 650 // Array of TStrings
 651   for(i=0;i<fIdN;i++) fIdSens[i] = gMC->VolId(fIdName[i]);
 652 //
 653 }
 654
 655 //_____________________________________________________________________________
 656 void AliITS::SetDefaults()
 657 {
 658   // sets the default segmentation, response, digit and raw cluster classes
 659
 660   printf("SetDefaults\n");
 661
 662   AliITSDetType *iDetType;
 663
 664
 665   //SPD
 666
 667   iDetType=DetType(0);
 668   if (!iDetType->GetSegmentationModel()) {
 669     AliITSsegmentationSPD *seg0=new AliITSsegmentationSPD(fITSgeom);
 670     SetSegmentationModel(0,seg0);
 671   }
 672   if (!iDetType->GetResponseModel()) {
 673      SetResponseModel(0,new AliITSresponseSPD());
 674   }
 675   // set digit and raw cluster classes to be used
 676
 677   const char *kData0=(iDetType->GetResponseModel())->DataType();
 678   if (strstr(kData0,"real")) {
 679       iDetType->ClassNames("AliITSdigit","AliITSRawClusterSPD");
 680   } else iDetType->ClassNames("AliITSdigitSPD","AliITSRawClusterSPD");
 681
 682   // SDD                                          //
 683   iDetType=DetType(1);
 684   if (!iDetType->GetResponseModel()) {
 685     SetResponseModel(1,new AliITSresponseSDD());
 686   }
 687   AliITSresponse *resp1=iDetType->GetResponseModel();
 688   if (!iDetType->GetSegmentationModel()) {
 689     AliITSsegmentationSDD *seg1=new AliITSsegmentationSDD(fITSgeom,resp1);
 690     SetSegmentationModel(1,seg1);
 691   }
 692   const char *kData1=(iDetType->GetResponseModel())->DataType();
 693   const char *kopt=iDetType->GetResponseModel()->ZeroSuppOption();
 694   if ((!strstr(kopt,"2D")) && (!strstr(kopt,"1D")) || strstr(kData1,"real") ) {
 695       iDetType->ClassNames("AliITSdigit","AliITSRawClusterSDD");
 696   } else iDetType->ClassNames("AliITSdigitSDD","AliITSRawClusterSDD");
 697
 698   // SSD
 699   iDetType=DetType(2);
 700   if (!iDetType->GetSegmentationModel()) {
 701     AliITSsegmentationSSD *seg2=new AliITSsegmentationSSD(fITSgeom);
 702     SetSegmentationModel(2,seg2);
 703   }
 704   if (!iDetType->GetResponseModel()) {
 705     SetResponseModel(2,new AliITSresponseSSD());
 706   }
 707   const char *kData2=(iDetType->GetResponseModel())->DataType();
 708   if (strstr(kData2,"real")) {
 709       iDetType->ClassNames("AliITSdigit","AliITSRawClusterSSD");
 710   } else iDetType->ClassNames("AliITSdigitSSD","AliITSRawClusterSSD");
 711
 712   if (kNTYPES>3) {
 713     Warning("SetDefaults","Only the three basic detector types are initialised!");
 714   }
 715
 716 }
 717
 718
 719
 720 //_____________________________________________________________________________
 721 void AliITS::SetDefaultSimulation()
 722 {
 723   // to be written
 724
 725 }
 726 //_____________________________________________________________________________
 727 void AliITS::SetDefaultClusterFinders()
 728 {
 729   // to be written
 730
 731 }
 732 //_____________________________________________________________________________
 733
 734 void AliITS::MakeTreeC(Option_t *option)
 735 {
 736   // create a separate tree to store the clusters
 737
 738 //  cout << "AliITS::MakeTreeC" << endl;
 739
 740      const char *optC = strstr(option,"C");
 741      if (optC && !fTreeC) fTreeC = new TTree("TC","Clusters in ITS");
 742      else return;
 743
 744      Int_t buffersize = 4000;
 745      char branchname[30];
 746
 747      const char *det[3] = {"SPD","SDD","SSD"};
 748
 749      char digclass[40];
 750      char clclass[40];
 751
 752      // one branch for Clusters per type of detector
 753      Int_t i;
 754      for (i=0; i<kNTYPES ;i++) {
 755        AliITSDetType *iDetType=DetType(i);
 756        iDetType->GetClassNames(digclass,clclass);
 757        // clusters
 758        (*fCtype)[i] = new TClonesArray(clclass,10000);
 759         if (kNTYPES==3) sprintf(branchname,"%sClusters%s",GetName(),det[i]);
 760         else  sprintf(branchname,"%sClusters%d",GetName(),i+1);
 761         if (fCtype   && fTreeC) {
 762            TreeC()->Branch(branchname,&((*fCtype)[i]), buffersize);
 763 //         cout << "Making Branch " << branchname;
 764 //         cout << " for Clusters of detector type " << i+1 << endl;
 765         }
 766      }
 767
 768 }
 769
 770 //_____________________________________________________________________________
 771 void AliITS::GetTreeC(Int_t event)
 772 {
 773
 774 //  cout << "AliITS::GetTreeC" << endl;
 775
 776   // get the clusters tree for this event and set the branch address
 777     char treeName[20];
 778     char branchname[30];
 779
 780     const char *det[3] = {"SPD","SDD","SSD"};
 781
 782     ResetClusters();
 783     if (fTreeC) {
 784           delete fTreeC;
 785     }
 786
 787     sprintf(treeName,"TreeC%d",event);
 788     fTreeC = (TTree*)gDirectory->Get(treeName);
 789
 790     TBranch *branch;
 791     if (fTreeC) {
 792         Int_t i;
 793         for (i=0; i<kNTYPES; i++) {
 794            if (kNTYPES==3) sprintf(branchname,"%sClusters%s",GetName(),det[i]);
 795            else  sprintf(branchname,"%sClusters%d",GetName(),i+1);
 796            if (fCtype) {
 797                 branch = fTreeC->GetBranch(branchname);
 798                 if (branch) branch->SetAddress(&((*fCtype)[i]));
 799             }
 800         }
 801     } else {
 802         Error("AliITS::GetTreeC",
 803                 "cannot find Clusters Tree for event:%d\n",event);
 804     }
 805
 806 }
 807 //_____________________________________________________________________________
 808 void AliITS::MakeBranch(Option_t* option, char *file)
 809 {
 810   //
 811   // Creates Tree branches for the ITS.
 812   //
 813   //
 814   Int_t buffersize = 4000;
 815   char branchname[30];
 816   sprintf(branchname,"%s",GetName());
 817
 818   AliDetector::MakeBranch(option,file);
 819
 820   const char *cD = strstr(option,"D");
 821   const char *cR = strstr(option,"R");
 822
 823   if (cD) {
 824   //
 825   // one branch for digits per type of detector
 826   //
 827    const char *det[3] = {"SPD","SDD","SSD"};
 828
 829    char digclass[40];
 830    char clclass[40];
 831
 832    Int_t i;
 833    for (i=0; i<kNTYPES ;i++) {
 834        AliITSDetType *iDetType=DetType(i);
 835        iDetType->GetClassNames(digclass,clclass);
 836        // digits
 837        if(!((*fDtype)[i])) (*fDtype)[i] = new TClonesArray(digclass,10000);
 838        else ResetDigits(i);
 839    }
 840
 841    for (i=0; i<kNTYPES ;i++) {
 842       if (kNTYPES==3) sprintf(branchname,"%sDigits%s",GetName(),det[i]);
 843       else  sprintf(branchname,"%sDigits%d",GetName(),i+1);
 844       if (fDtype && gAlice->TreeD()) {
 845         gAlice->MakeBranchInTree(gAlice->TreeD(),
 846                          branchname, &((*fDtype)[i]), buffersize, file);
 847 //          cout << "Making Branch " << branchname;
 848 //          cout << " for digits of type "<< i+1 << endl;
 849       }
 850     }
 851   }
 852
 853   if (cR) {
 854   //
 855   // only one branch for rec points for all detector types
 856   //
 857     sprintf(branchname,"%sRecPoints",GetName());
 858
 859     if(!fRecPoints) fRecPoints=new TClonesArray("AliITSRecPoint",10000);
 860
 861     if (fRecPoints && gAlice->TreeR()) {
 862       gAlice->MakeBranchInTree(gAlice->TreeR(),
 863                                branchname, &fRecPoints, buffersize, file) ;
 864 //      cout << "Making Branch " << branchname;
 865 //      cout << " for reconstructed space points" << endl;
 866     }
 867   }
 868 }
 869
 870 //___________________________________________
 871 void AliITS::SetTreeAddress()
 872 {
 873
 874   // Set branch address for the Trees.
 875
 876   char branchname[30];
 877   AliDetector::SetTreeAddress();
 878
 879   const char *det[3] = {"SPD","SDD","SSD"};
 880
 881   TBranch *branch;
 882   TTree *treeD = gAlice->TreeD();
 883   TTree *treeR = gAlice->TreeR();
 884
 885   Int_t i;
 886   if (treeD) {
 887       for (i=0; i<kNTYPES; i++) {
 888           if (kNTYPES==3) sprintf(branchname,"%sDigits%s",GetName(),det[i]);
 889           else  sprintf(branchname,"%sDigits%d",GetName(),i+1);
 890           if (fDtype) {
 891               branch = treeD->GetBranch(branchname);
 892               if (branch) branch->SetAddress(&((*fDtype)[i]));
 893           }
 894       }
 895   }
 896
 897
 898   if (treeR) {
 899     sprintf(branchname,"%sRecPoints",GetName());
 900       branch = treeR->GetBranch(branchname);
 901       if (branch) branch->SetAddress(&fRecPoints);
 902   }
 903
 904
 905 }
 906
 907 //____________________________________________________________________________
 908 void AliITS::InitModules(Int_t size,Int_t &nmodules){
 909
 910   //initialize the modules array
 911
 912   if(fITSmodules){
 913       fITSmodules->Delete();
 914       delete fITSmodules;
 915   }
 916
 917     Int_t nl,indexMAX,index;
 918
 919     if(size<=0){ // default to using data stored in AliITSgeom
 920         if(fITSgeom==0) {
 921             Error("AliITS::InitModules",
 922                 "in AliITS::InitModule fITSgeom not defined\n");
 923             return;
 924         } // end if fITSgeom==0
 925         nl = fITSgeom->GetNlayers();
 926         indexMAX = fITSgeom->GetModuleIndex(nl,fITSgeom->GetNladders(nl),
 927                                             fITSgeom->GetNdetectors(nl))+1;
 928         nmodules = indexMAX;
 929         fITSmodules = new TObjArray(indexMAX);
 930         for(index=0;index<indexMAX;index++){
 931                 fITSmodules->AddAt( new AliITSmodule(index),index);
 932         } // end for index
 933     }else{
 934         fITSmodules = new TObjArray(size);
 935         for(index=0;index<size;index++) {
 936             fITSmodules->AddAt( new AliITSmodule(index),index);
 937         }
 938
 939         nmodules = size;
 940     } // end i size<=0
 941 }
 942
 943 //____________________________________________________________________________
 944 void AliITS::FillModules(Int_t evnt,Int_t bgrev,Int_t nmodules,Option_t *option,Text_t *filename){
 945
 946   // fill the modules with the sorted by module hits; add hits from background
 947   // if option=Add
 948
 949
 950     static TTree *trH1;                 //Tree with background hits
 951     static TClonesArray *fHits2;        //List of hits for one track only
 952
 953     static Bool_t first=kTRUE;
 954     static TFile *file;
 955     const char *addBgr = strstr(option,"Add");
 956
 957
 958     if (addBgr ) {
 959         if(first) {
 960 //          cout<<"filename "<<filename<<endl;
 961             file=new TFile(filename);
 962 //          cout<<"I have opened "<<filename<<" file "<<endl;
 963             fHits2     = new TClonesArray("AliITShit",1000  );
 964         }
 965         first=kFALSE;
 966         file->cd();
 967         file->ls();
 968         // Get Hits Tree header from file
 969         if(fHits2) fHits2->Clear();
 970         if(trH1) delete trH1;
 971         trH1=0;
 972
 973         char treeName[20];
 974         sprintf(treeName,"TreeH%d",bgrev);
 975         trH1 = (TTree*)gDirectory->Get(treeName);
 976         //printf("TrH1 %p of treename %s for event %d \n",trH1,treeName,bgrev);
 977
 978         if (!trH1) {
 979             Error("AliITS::FillModules",
 980             "cannot find Hits Tree for event:%d\n",bgrev);
 981         }
 982         // Set branch addresses
 983         TBranch *branch;
 984         char branchname[20];
 985         sprintf(branchname,"%s",GetName());
 986         if (trH1 && fHits2) {
 987             branch = trH1->GetBranch(branchname);
 988             if (branch) branch->SetAddress(&fHits2);
 989         }
 990
 991         // test
 992         //Int_t ntracks1 =(Int_t)TrH1->GetEntries();
 993         //printf("background - ntracks1 - %d\n",ntracks1);
 994    }
 995
 996     //Int_t npart = gAlice->GetEvent(evnt);
 997     //if(npart<=0) return;
 998     TClonesArray *itsHits = this->Hits();
 999     Int_t lay,lad,det,index;
1000     AliITShit *itsHit=0;
1001     AliITSmodule *mod=0;
1002
1003     TTree *iTH = gAlice->TreeH();
1004     Int_t ntracks =(Int_t) iTH->GetEntries();
1005
1006     Int_t t,h;
1007     for(t=0; t<ntracks; t++){
1008         gAlice->ResetHits();
1009         iTH->GetEvent(t);
1010         Int_t nhits = itsHits->GetEntriesFast();
1011         //printf("nhits %d\n",nhits);
1012         if (!nhits) continue;
1013         for(h=0; h<nhits; h++){
1014             itsHit = (AliITShit *)itsHits->UncheckedAt(h);
1015             itsHit->GetDetectorID(lay,lad,det);
1016             // temporarily index=det-1 !!!
1017             if(fITSgeom) index = fITSgeom->GetModuleIndex(lay,lad,det);
1018             else index=det-1;
1019             //
1020             mod = this->GetModule(index);
1021             mod->AddHit(itsHit,t,h);
1022         } // end loop over hits
1023     } // end loop over tracks
1024
1025     // open the file with background
1026
1027     if (addBgr ) {
1028           Int_t track,i;
1029           ntracks =(Int_t)trH1->GetEntries();
1030             //printf("background - ntracks1 %d\n",ntracks);
1031             //printf("background - Start loop over tracks \n");
1032             //   Loop over tracks
1033
1034             for (track=0; track<ntracks; track++) {
1035
1036                 if (fHits2)       fHits2->Clear();
1037                 trH1->GetEvent(track);
1038                 //   Loop over hits
1039                 for(i=0;i<fHits2->GetEntriesFast();++i) {
1040
1041                     itsHit=(AliITShit*) (*fHits2)[i];
1042                     itsHit->GetDetectorID(lay,lad,det);
1043                     // temporarily index=det-1 !!!
1044                     if(fITSgeom) index = fITSgeom->GetModuleIndex(lay,lad,det);
1045                     else index=det-1;
1046                     //
1047                     mod = this->GetModule(index);
1048                     mod->AddHit(itsHit,track,i);
1049                }  // end loop over hits
1050             } // end loop over tracks
1051
1052             TTree *fAli=gAlice->TreeK();
1053             TFile *fileAli=0;
1054
1055             if (fAli) fileAli =fAli->GetCurrentFile();
1056             fileAli->cd();
1057
1058     } // end if add
1059
1060     //gObjectTable->Print();
1061
1062 }
1063
1064 //____________________________________________________________________________
1065
1066 void AliITS::SDigits2Digits()
1067 {
1068
1069   AliITSgeom *geom = GetITSgeom();
1070
1071   // SPD
1072   AliITSDetType *iDetType;
1073   iDetType=DetType(0);
1074   AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
1075   AliITSresponseSPD *res0 = (AliITSresponseSPD*)iDetType->GetResponseModel();
1076   AliITSsimulationSPD *sim0=new AliITSsimulationSPD(seg0,res0);
1077   SetSimulationModel(0,sim0);
1078   // test
1079   // printf("SPD dimensions %f %f \n",seg0->Dx(),seg0->Dz());
1080   // printf("SPD npixels %d %d \n",seg0->Npz(),seg0->Npx());
1081   // printf("SPD pitches %d %d \n",seg0->Dpz(0),seg0->Dpx(0));
1082   // end test
1083   //
1084   // SDD
1085   //Set response functions
1086   // SDD compression param: 2 fDecrease, 2fTmin, 2fTmax or disable, 2 fTolerance
1087
1088   iDetType=DetType(1);
1089   AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
1090   if (!res1) {
1091     res1=new AliITSresponseSDD();
1092     SetResponseModel(1,res1);
1093   }
1094   Float_t noise, baseline;
1095   res1->GetNoiseParam(noise,baseline);
1096   Float_t noise_after_el =  res1->GetNoiseAfterElectronics();
1097   Float_t fCutAmp = baseline + 2.*noise_after_el;
1098   Int_t cp[8]={0,0,(int)fCutAmp,(int)fCutAmp,0,0,0,0}; //1D
1099   res1->SetCompressParam(cp);
1100   AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
1101   if (!seg1) {
1102     seg1 = new AliITSsegmentationSDD(geom,res1);
1103     SetSegmentationModel(1,seg1);
1104   }
1105   AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1);
1106   SetSimulationModel(1,sim1);
1107
1108   // SSD
1109   iDetType=DetType(2);
1110   AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
1111   AliITSresponseSSD *res2 = (AliITSresponseSSD*)iDetType->GetResponseModel();
1112   res2->SetSigmaSpread(3.,2.);
1113   AliITSsimulationSSD *sim2=new AliITSsimulationSSD(seg2,res2);
1114   SetSimulationModel(2,sim2);
1115
1116   cerr<<"Digitizing ITS...\n";
1117
1118   TStopwatch timer;
1119   timer.Start();
1120   HitsToDigits(0,0,-1," ","All"," ");
1121   timer.Stop(); timer.Print();
1122
1123   delete sim0;
1124   delete sim1;
1125   delete sim2;
1126 }
1127
1128
1129 //____________________________________________________________________________
1130 void AliITS::HitsToDigits(Int_t evNumber,Int_t bgrev,Int_t size, Option_t *option, Option_t *opt,Text_t *filename)
1131 {
1132     // keep galice.root for signal and name differently the file for
1133     // background when add! otherwise the track info for signal will be lost !
1134
1135    // the condition below will disappear when the geom class will be
1136    // initialised for all versions - for the moment it is only for v5 !
1137    // 7 is the SDD beam test version
1138    Int_t ver = this->IsVersion();
1139    if(ver!=5 && ver!=7 && ver!=8 && ver!=9) return;
1140
1141    const char *all = strstr(opt,"All");
1142    const char *det[3] = {strstr(opt,"SPD"),strstr(opt,"SDD"),strstr(opt,"SSD")};
1143 //   cout<<" 1 AliITS "<<endl;
1144    Int_t nmodules;
1145    InitModules(size,nmodules);
1146 //   cout<<" 2 AliITS "<<endl;
1147    FillModules(evNumber,bgrev,nmodules,option,filename);
1148 //   cout<<" 3 AliITS "<<endl;
1149
1150    //TBranch *branch;
1151    AliITSsimulation* sim;
1152    //TObjArray *branches=gAlice->TreeD()->GetListOfBranches();
1153    AliITSgeom *geom = GetITSgeom();
1154
1155    Int_t id,module;
1156 //   Int_t lay, lad, detect;
1157    Int_t first,last;
1158    for (id=0;id<kNTYPES;id++) {
1159         if (!all && !det[id]) continue;
1160         //branch = (TBranch*)branches->UncheckedAt(id);
1161         AliITSDetType *iDetType=DetType(id);
1162         sim = (AliITSsimulation*)iDetType->GetSimulationModel();
1163
1164         if (!sim) {
1165            Error("HitsToDigits","The simulation class was not instantiated!");
1166            exit(1);
1167            // or SetDefaultSimulation();
1168         }
1169
1170         if(geom) {
1171           first = geom->GetStartDet(id);
1172           last = geom->GetLastDet(id);
1173         } else first=last=0;
1174 //      cout << "det type " << id << " first, last "<< first << last << endl;
1175         for(module=first;module<=last;module++) {
1176             AliITSmodule *mod = (AliITSmodule *)fITSmodules->At(module);
1177 /*
1178             geom->GetModuleId(module,lay, lad, detect);
1179             if ( lay == 6 )
1180               ((AliITSsegmentationSSD*)(((AliITSsimulationSSD*)sim)->GetSegmentation()))->SetLayer(6);
1181             if ( lay == 5 )
1182               ((AliITSsegmentationSSD*)(((AliITSsimulationSSD*)sim)->GetSegmentation()))->SetLayer(5);
1183 */
1184             sim->DigitiseModule(mod,module,evNumber);
1185             // fills all branches - wasted disk space
1186             gAlice->TreeD()->Fill();
1187             ResetDigits();
1188             // try and fill only the branch
1189             //branch->Fill();
1190             //ResetDigits(id);
1191         } // loop over modules
1192    } // loop over detector types
1193
1194    ClearModules();
1195
1196 //   Int_t nentries=(Int_t)
1197    gAlice->TreeD()->GetEntries();
1198 //   cout << "nentries in TreeD" << nentries << endl;
1199
1200    char hname[30];
1201    sprintf(hname,"TreeD%d",evNumber);
1202    gAlice->TreeD()->Write(hname,TObject::kOverwrite);
1203    // reset tree
1204    gAlice->TreeD()->Reset();
1205
1206 }
1207
1208
1209 //____________________________________________________________________________
1210 void AliITS::DigitsToRecPoints(Int_t evNumber,Int_t lastentry,Option_t *opt)
1211 {
1212   // cluster finding and reconstruction of space points
1213
1214    // the condition below will disappear when the geom class will be
1215    // initialised for all versions - for the moment it is only for v5 !
1216    // 7 is the SDD beam test version
1217    Int_t ver = this->IsVersion();
1218    if(ver!=5 && ver!=8 && ver!=9) return;
1219
1220    const char *all = strstr(opt,"All");
1221    const char *det[3] = {strstr(opt,"SPD"),strstr(opt,"SDD"),strstr(opt,"SSD")};
1222
1223    static Bool_t first=kTRUE;
1224    if (!TreeC() && first) {
1225        MakeTreeC("C");
1226        first=kFALSE;
1227    }
1228
1229    TTree *treeC=TreeC();
1230
1231
1232    //TBranch *branch;
1233    AliITSClusterFinder* rec;
1234
1235    //TObjArray *branches=gAlice->TreeR()->GetListOfBranches();
1236    AliITSgeom *geom = GetITSgeom();
1237
1238    Int_t id,module;
1239    for (id=0;id<kNTYPES;id++) {
1240         if (!all && !det[id]) continue;
1241         //branch = (TBranch*)branches->UncheckedAt(id);
1242         AliITSDetType *iDetType=DetType(id);
1243         rec = (AliITSClusterFinder*)iDetType->GetReconstructionModel();
1244         if (!rec) {
1245            Error("DigitsToRecPoints","The cluster finder class was not instantiated!");
1246            exit(1);
1247            // or SetDefaultClusterFinders();
1248         }
1249         TClonesArray *itsDigits  = this->DigitsAddress(id);
1250
1251         Int_t first,last;
1252         if(geom) {
1253           first = geom->GetStartDet(id);
1254           last = geom->GetLastDet(id);
1255         } else first=last=0;
1256         printf("first module - last module %d %d\n",first,last);
1257         for(module=first;module<=last;module++) {
1258               this->ResetDigits();
1259               if (all) gAlice->TreeD()->GetEvent(lastentry+module);
1260               else gAlice->TreeD()->GetEvent(lastentry+(module-first));
1261               Int_t ndigits = itsDigits->GetEntriesFast();
1262               if (ndigits) rec->FindRawClusters(module);
1263               gAlice->TreeR()->Fill();
1264               ResetRecPoints();
1265               treeC->Fill();
1266               ResetClusters();
1267               // try and fill only the branch
1268               //branch->Fill();
1269               //ResetRecPoints(id);
1270         } // loop over modules
1271    } // loop over detector types
1272
1273
1274 //   Int_t nentries=(Int_t)
1275    gAlice->TreeR()->GetEntries();
1276 //   Int_t ncentries=(Int_t)
1277    treeC->GetEntries();
1278 //   cout << " nentries ncentries " << nentries << ncentries <<  endl;
1279
1280    char hname[30];
1281    sprintf(hname,"TreeR%d",evNumber);
1282    gAlice->TreeR()->Write(hname,TObject::kOverwrite);
1283    // reset tree
1284    gAlice->TreeR()->Reset();
1285
1286    sprintf(hname,"TreeC%d",evNumber);
1287    treeC->Write(hname,TObject::kOverwrite);
1288    treeC->Reset();
1289 }
1290
1291
1292 //____________________________________________________________________________
1293 void AliITS::HitsToFastRecPoints(Int_t evNumber,Int_t bgrev,Int_t size,
1294 Option_t *option,Option_t *opt,Text_t *filename)
1295 {
1296     // keep galice.root for signal and name differently the file for
1297     // background when add! otherwise the track info for signal will be lost !
1298
1299
1300    // the condition below will disappear when the geom class will be
1301    // initialised for all versions - for the moment it is only for v5 !
1302    Int_t ver = this->IsVersion();
1303    if(ver!=5 && ver!=8 && ver!=9) return;
1304    //if(ver!=5) return;
1305
1306    const char *all = strstr(opt,"All");
1307    const char *det[3] = {strstr(opt,"SPD"),strstr(opt,"SDD"),strstr(opt,"SSD")};
1308
1309    Int_t nmodules;
1310    InitModules(size,nmodules);
1311    FillModules(evNumber,bgrev,nmodules,option,filename);
1312
1313
1314    AliITSsimulation* sim;
1315    AliITSgeom *geom = GetITSgeom();
1316
1317    TRandom *random=new TRandom[9];
1318    random[0].SetSeed(111);
1319    random[1].SetSeed(222);
1320    random[2].SetSeed(333);
1321    random[3].SetSeed(444);
1322    random[4].SetSeed(555);
1323    random[5].SetSeed(666);
1324    random[6].SetSeed(777);
1325    random[7].SetSeed(888);
1326    random[8].SetSeed(999);
1327
1328
1329    Int_t id,module;
1330    for (id=0;id<kNTYPES;id++) {
1331         if (!all && !det[id]) continue;
1332         AliITSDetType *iDetType=DetType(id);
1333         sim = (AliITSsimulation*)iDetType->GetSimulationModel();
1334         if (!sim) {
1335            Error("HitsToFastPoints","The simulation class was not instantiated!");
1336            exit(1);
1337            // or SetDefaultSimulation();
1338         }
1339
1340         Int_t first,last;
1341         if(geom) {
1342           first = geom->GetStartDet(id);
1343           last = geom->GetLastDet(id);
1344         } else first=last=0;
1345         printf("first module - last module %d %d\n",first,last);
1346         for(module=first;module<=last;module++) {
1347             AliITSmodule *mod = (AliITSmodule *)fITSmodules->At(module);
1348             sim->CreateFastRecPoints(mod,module,random);
1349             gAlice->TreeR()->Fill();
1350             ResetRecPoints();
1351         } // loop over modules
1352    } // loop over detector types
1353
1354
1355    ClearModules();
1356
1357    //Int_t nentries=(Int_t)gAlice->TreeR()->GetEntries();
1358
1359    char hname[30];
1360    sprintf(hname,"TreeR%d",evNumber);
1361    gAlice->TreeR()->Write(hname,TObject::kOverwrite);
1362    // reset tree
1363    gAlice->TreeR()->Reset();
1364
1365    delete [] random;
1366
1367 }
1368
1369 //________________________________________________________________
1370 AliITStrack  AliITS::Tracking(AliITStrack &track, AliITStrack *reference,TObjArray *fastpoints, Int_t
1371 **vettid, Bool_t flagvert,  AliITSRad *rl  ) {
1372
1373 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
1374
1375
1376   TList *list= new TList();
1377
1378   AliITStrack tr(track);
1379
1380   list->AddLast(&tr);
1381
1382   Double_t Pt=(tr).GetPt();
1383 //  cout << "\n Pt = " << Pt <<"\n";  //stampa
1384
1385   AliITStracking obj(list, reference, this, fastpoints,TMath::Abs(Pt),vettid, flagvert, rl);
1386   list->Delete();
1387   delete list;
1388
1389   Int_t itot=-1;
1390   TVector VecTotLabref(18);
1391   Int_t lay, k;
1392   for(lay=5; lay>=0; lay--) {
1393     TVector VecLabref(3);
1394     VecLabref=(*reference).GetLabTrack(lay);
1395     Float_t ClustZ=(*reference).GetZclusterTrack( lay);   //aggiunta il 5-3-2001
1396     for(k=0; k<3; k++){ //{itot++; VecTotLabref(itot)=VecLabref(k);} //cambiata il 5-3-2001
1397                 Int_t lpp=(Int_t)VecLabref(k);
1398                 if(lpp>=0) {
1399                   TParticle *p=(TParticle*) gAlice->Particle(lpp);
1400                   Int_t pcode=p->GetPdgCode();
1401                   if(pcode==11) VecLabref(k)=p->GetFirstMother();
1402                 }
1403     itot++; VecTotLabref(itot)=VecLabref(k);
1404     if(VecLabref(k)==0. && ClustZ == 0.) VecTotLabref(itot) =-3.; }
1405   }
1406   Long_t labref;
1407   Int_t freq;
1408   (*reference).Search(VecTotLabref, labref, freq);
1409
1410   //if(freq < 4) labref=-labref;
1411   //if(freq < 6) labref=-labref;        // cinque - sei
1412   if(freq < 5) labref=-labref;        // cinque - sei
1413   (*reference).SetLabel(labref);
1414
1415   return *reference;
1416
1417 }
1418
1419
1420
1421 //________________________________________________________________
1422
1423
1424
1425 void AliITS::DoTracking(Int_t evNumber, Int_t min_t, Int_t max_t, TFile *file, Bool_t flagvert) {
1426
1427 //   ex macro for tracking ITS
1428
1429   printf("begin DoTracking - file %p\n",file);
1430
1431   //const char *pname="75x40_100x60";
1432
1433  Int_t imax=200,jmax=450;
1434   AliITSRad *rl = new AliITSRad(imax,jmax);
1435   //cout<<" dopo costruttore AliITSRad\n"; getchar();
1436
1437   struct GoodTrack {
1438     Int_t lab,code;
1439     Float_t px,py,pz,x,y,z,pxg,pyg,pzg,ptg;
1440     Bool_t flag;
1441   };
1442
1443
1444   gAlice->GetEvent(0);
1445
1446   //  AliKalmanTrack *kkprov;
1447   //  kkprov->SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
1448
1449   /*  //modificato il 26-4-2001
1450   AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
1451   AliTPCParam *digp = (AliTPCParam*)file->Get(pname);
1452   if (digp!=0) TPC->SetParam(digp);
1453   */
1454   TFile *cf=TFile::Open("AliTPCclusters.root");
1455   AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
1456   if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();}
1457
1458
1459   GoodTrack gt[15000];
1460   Int_t ngood=0;
1461   ifstream in("itsgood_tracks");
1462
1463   cerr<<"Reading itsgood tracks...\n";
1464   while (in>>gt[ngood].lab>>gt[ngood].code
1465           >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
1466           >>gt[ngood].x  >>gt[ngood].y >>gt[ngood].z
1467           >>gt[ngood].pxg  >>gt[ngood].pyg >>gt[ngood].pzg
1468           >>gt[ngood].ptg >>gt[ngood].flag) {
1469     ngood++;
1470     cerr<<ngood<<'\r';
1471     if (ngood==15000) {
1472       cerr<<"Too many good tracks !\n";
1473       break;
1474     }
1475   }
1476   if (!in.eof()) cerr<<"Read error (itsgood_tracks) !\n";
1477
1478
1479 // Load tracks
1480  // TFile *tf=TFile::Open("tpctracks.root"); //commentato 26-4-2001
1481   TFile *tf=TFile::Open("AliTPCtracks.root");
1482   if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return ;}
1483   TObjArray tracks(200000);
1484   TTree *tracktree=(TTree*)tf->Get("TreeT");
1485   TBranch *tbranch=tracktree->GetBranch("tracks");
1486   Int_t nentr=(Int_t)tracktree->GetEntries();
1487   Int_t kk;
1488   /*  commentato il 26-4-2001
1489   for (kk=0; kk<nentr; kk++) {
1490     AliTPCtrack *iotrack=new AliTPCtrack;
1491     tbranch->SetAddress(&iotrack);
1492     tracktree->GetEvent(kk);
1493     tracks.AddLast(iotrack);
1494   }
1495   */
1496    AliTPCtrack *iotracktpc=0;
1497   for (kk=0; kk<nentr; kk++) {
1498     iotracktpc=new AliTPCtrack;
1499     tbranch->SetAddress(&iotracktpc);
1500     tracktree->GetEvent(kk);
1501     tracks.AddLast(iotracktpc);
1502   }
1503   tf->Close();
1504
1505
1506   Int_t nt = tracks.GetEntriesFast();
1507   cerr<<"Number of found tracks "<<nt<<endl;
1508
1509   TVector DataOut(9);
1510   Int_t kkk=0;
1511
1512   Double_t ptg=0.,pxg=0.,pyg=0.,pzg=0.;
1513
1514   //////////////////////////////  good tracks definition in TPC  ////////////////////////////////
1515
1516   ofstream out1 ("AliITSTrag.out");
1517   Int_t i;
1518   for (i=0; i<ngood; i++) out1 << gt[i].ptg << "\n";
1519   out1.close();
1520
1521
1522   TVector vec(5);
1523   TTree *TR=gAlice->TreeR();
1524   Int_t nent=(Int_t)TR->GetEntries();
1525   TClonesArray  *recPoints = RecPoints();
1526   Int_t numbpoints;
1527   Int_t totalpoints=0;
1528   Int_t *np = new Int_t[nent];
1529   Int_t **vettid = new Int_t* [nent];
1530   Int_t mod;
1531   for (mod=0; mod<nent; mod++) {
1532     vettid[mod]=0;
1533     this->ResetRecPoints();
1534     //gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
1535     gAlice->TreeR()->GetEvent(mod); //first entry in TreeR is empty
1536     numbpoints = recPoints->GetEntries();
1537     totalpoints+=numbpoints;
1538     np[mod] = numbpoints;
1539   //cout<<" mod = "<<mod<<"   numbpoints = "<<numbpoints<<"\n"; getchar();
1540     vettid[mod] = new Int_t[numbpoints];
1541     Int_t ii;
1542     for (ii=0;ii<numbpoints; ii++) *(vettid[mod]+ii)=0;
1543   }
1544
1545   AliTPCtrack *track=0;
1546
1547
1548   if(min_t < 0) {min_t = 0; max_t = nt-1;}
1549
1550 /*
1551   ///////////////////////////////// Definition of vertex end its error ////////////////////////////
1552   ////////////////////////// In the future it will be given by a method ///////////////////////////
1553   Double_t Vx=0.;
1554   Double_t Vy=0.;
1555   Double_t Vz=0.;
1556
1557   Float_t sigmavx=0.0050;      // 50  microns
1558   Float_t sigmavy=0.0050;      // 50  microns
1559   Float_t sigmavz=0.010;       // 100 microns
1560
1561   //Vx+=gRandom->Gaus(0,sigmavx);  Vy+=gRandom->Gaus(0,sigmavy);  Vz+=gRandom->Gaus(0,sigmavz);
1562   TVector vertex(3), ervertex(3)
1563   vertex(0)=Vx; vertex(1)=Vy; vertex(2)=Vz;
1564   ervertex(0)=sigmavx;  ervertex(1)=sigmavy;  ervertex(2)=sigmavz;
1565   /////////////////////////////////////////////////////////////////////////////////////////////////
1566 */
1567
1568
1569   TTree tracktree1("TreeT","Tree with ITS tracks");
1570   AliITSiotrack *iotrack=0;
1571   tracktree1.Branch("ITStracks","AliITSiotrack",&iotrack,32000,0);
1572
1573   ofstream out ("AliITSTra.out");
1574   //ofstream outprova ("AliITSprova.out"); //commentato il 26-4-2001
1575
1576   Int_t j;
1577   for (j=min_t; j<=max_t; j++) {
1578     track=(AliTPCtrack*)tracks.UncheckedAt(j);
1579     Int_t flaglab=0;
1580     if (!track) continue;
1581     ////// elimination of not good tracks ////////////
1582     Int_t ilab=TMath::Abs(track->GetLabel());
1583     Int_t iii;
1584     for (iii=0;iii<ngood;iii++) {
1585          //cout<<" ilab, gt[iii].lab = "<<ilab<<" "<<gt[iii].lab<<"\n"; getchar();
1586       if (ilab==gt[iii].lab) {
1587         flaglab=1;
1588         ptg=gt[iii].ptg;
1589         pxg=gt[iii].pxg;
1590         pyg=gt[iii].pyg;
1591         pzg=gt[iii].pzg;
1592         break;
1593       }
1594     }
1595          //cout<<" j flaglab =  " <<j<<" "<<flaglab<<"\n";  getchar();
1596     if (!flaglab) continue;
1597          //cout<<" j =  " <<j<<"\n";  getchar();
1598   /*
1599     ////// old propagation to the end of TPC //////////////
1600     Double_t xk=76.;
1601     track->PropagateTo(xk);
1602     xk-=0.11;
1603     track->PropagateTo(xk,42.7,2.27); //C
1604     xk-=2.6;
1605     track->PropagateTo(xk,36.2,1.98e-3); //C02
1606     xk-=0.051;
1607     track->PropagateTo(xk,42.7,2.27); //C
1608     ///////////////////////////////////////////////////
1609          */
1610
1611          ////// new propagation to the end of TPC //////////////
1612     Double_t xk=77.415;
1613     track->PropagateTo(xk, 28.94, 1.204e-3);     //Ne
1614          xk -=0.01;
1615     track->PropagateTo(xk, 44.77, 1.71);         //Tedlar
1616          xk -=0.04;
1617     track->PropagateTo(xk, 44.86, 1.45);         //Kevlar
1618          xk -=2.0;
1619     track->PropagateTo(xk, 41.28, 0.029);        //Nomex
1620     xk-=16;
1621     track->PropagateTo(xk,36.2,1.98e-3); //C02
1622          xk -=0.01;
1623     track->PropagateTo(xk, 24.01, 2.7);  //Al
1624          xk -=0.01;
1625     track->PropagateTo(xk, 44.77, 1.71);         //Tedlar
1626          xk -=0.04;
1627     track->PropagateTo(xk, 44.86, 1.45);         //Kevlar
1628          xk -=0.5;
1629     track->PropagateTo(xk, 41.28, 0.029);        //Nomex
1630
1631        ///////////////////////////////////////////////////////////////
1632
1633    ///////////////////////////////////////////////////////////////
1634     AliITStrack trackITS(*track);
1635     AliITStrack result(*track);
1636     AliITStrack primarytrack(*track);
1637
1638 ///////////////////////////////////////////////////////////////////////////////////////////////
1639          TVector Vgeant(3);
1640          Vgeant=result.GetVertex();
1641
1642   // Definition of Dv and Zv for vertex constraint
1643      Double_t sigmaDv=0.0050;  Double_t sigmaZv=0.010;
1644     //Double_t sigmaDv=0.0015;  Double_t sigmaZv=0.0015;
1645         Double_t uniform= gRandom->Uniform();
1646         Double_t signdv;
1647         if(uniform<=0.5) signdv=-1.;
1648            else
1649                  signdv=1.;
1650
1651         Double_t Vr=TMath::Sqrt(Vgeant(0)*Vgeant(0)+ Vgeant(1)*Vgeant(1));
1652           Double_t Dv=gRandom->Gaus(signdv*Vr,(Float_t)sigmaDv);
1653     Double_t Zv=gRandom->Gaus(Vgeant(2),(Float_t)sigmaZv);
1654
1655   //cout<<" Dv e Zv = "<<Dv<<" "<<Zv<<"\n";
1656     trackITS.SetDv(Dv);  trackITS.SetZv(Zv);
1657     trackITS.SetsigmaDv(sigmaDv); trackITS.SetsigmaZv(sigmaZv);
1658     result.SetDv(Dv);  result.SetZv(Zv);
1659     result.SetsigmaDv(sigmaDv); result.SetsigmaZv(sigmaZv);
1660     primarytrack.SetDv(Dv);  primarytrack.SetZv(Zv);
1661     primarytrack.SetsigmaDv(sigmaDv); primarytrack.SetsigmaZv(sigmaZv);
1662
1663 /////////////////////////////////////////////////////////////////////////////////////////////////
1664
1665     primarytrack.PrimaryTrack(rl);
1666     TVector  d2=primarytrack.Getd2();
1667     TVector  tgl2=primarytrack.Gettgl2();
1668     TVector  dtgl=primarytrack.Getdtgl();
1669     trackITS.Setd2(d2); trackITS.Settgl2(tgl2);  trackITS.Setdtgl(dtgl);
1670     result.Setd2(d2); result.Settgl2(tgl2);  result.Setdtgl(dtgl);
1671          /*
1672     trackITS.SetVertex(vertex); trackITS.SetErrorVertex(ervertex);
1673     result.SetVertex(vertex);   result.SetErrorVertex(ervertex);
1674     */
1675
1676     Tracking(trackITS,&result,recPoints,vettid, flagvert,rl);
1677
1678     // cout<<" progressive track number = "<<j<<"\r";
1679    // cout<<j<<"\r";
1680     Int_t NumofCluster=result.GetNumClust();
1681 //    cout<<" progressive track number = "<<j<<"\n";    // stampa
1682     Long_t labITS=result.GetLabel();
1683 //    cout << " ITS track label = " << labITS << "\n";  // stampa
1684     int lab=track->GetLabel();
1685 //    cout << " TPC track label = " << lab <<"\n";      // stampa
1686
1687
1688 //propagation to vertex
1689
1690     Double_t rbeam=3.;
1691
1692     result.Propagation(rbeam);
1693
1694          Double_t C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44;
1695          result.GetCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);
1696
1697     Double_t pt=TMath::Abs(result.GetPt());
1698     Double_t Dr=result.GetD();
1699     Double_t Z=result.GetZ();
1700     Double_t tgl=result.GetTgl();
1701     Double_t C=result.GetC();
1702     Double_t Cy=C/2.;
1703     Double_t Dz=Z-(tgl/Cy)*TMath::ASin(result.arga(rbeam));
1704          Dz-=Vgeant(2);
1705
1706          // cout<<" Dr e dz alla fine = "<<Dr<<" "<<Dz<<"\n"; getchar();
1707     Double_t phi=result.Getphi();
1708     Double_t phivertex = phi - TMath::ASin(result.argA(rbeam));
1709     Double_t duepi=2.*TMath::Pi();
1710     if(phivertex>duepi) phivertex-=duepi;
1711     if(phivertex<0.) phivertex+=duepi;
1712     Double_t Dtot=TMath::Sqrt(Dr*Dr+Dz*Dz);
1713
1714 //////////////////////////////////////////////////////////////////////////////////////////
1715
1716     Int_t idmodule,idpoint;
1717          if(NumofCluster >=5)  {            // cinque - sei
1718          //if(NumofCluster ==6)  {            // cinque - sei
1719
1720
1721       AliITSiotrack outtrack;
1722
1723       iotrack=&outtrack;
1724
1725       iotrack->SetStatePhi(phi);
1726       iotrack->SetStateZ(Z);
1727       iotrack->SetStateD(Dr);
1728       iotrack->SetStateTgl(tgl);
1729       iotrack->SetStateC(C);
1730                 Double_t radius=result.Getrtrack();
1731                 iotrack->SetRadius(radius);
1732                 Int_t charge;
1733                 if(C>0.) charge=-1;  else charge=1;
1734                 iotrack->SetCharge(charge);
1735
1736
1737
1738       iotrack->SetCovMatrix(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);
1739
1740       Double_t px=pt*TMath::Cos(phivertex);
1741       Double_t py=pt*TMath::Sin(phivertex);
1742       Double_t pz=pt*tgl;
1743
1744       Double_t xtrack=Dr*TMath::Sin(phivertex);
1745       Double_t ytrack=Dr*TMath::Cos(phivertex);
1746       Double_t ztrack=Dz+Vgeant(2);
1747
1748
1749       iotrack->SetPx(px);
1750       iotrack->SetPy(py);
1751       iotrack->SetPz(pz);
1752       iotrack->SetX(xtrack);
1753       iotrack->SetY(ytrack);
1754       iotrack->SetZ(ztrack);
1755       iotrack->SetLabel(labITS);
1756
1757       Int_t il;
1758                 for(il=0;il<6; il++){
1759                   iotrack->SetIdPoint(il,result.GetIdPoint(il));
1760                   iotrack->SetIdModule(il,result.GetIdModule(il));
1761                 }
1762       tracktree1.Fill();
1763
1764    //cout<<" labITS = "<<labITS<<"\n";
1765         //cout<<" phi z Dr tgl C = "<<phi<<" "<<Z<<" "<<Dr<<" "<<tgl<<" "<<C<<"\n";  getchar();
1766
1767      DataOut(kkk) = ptg; kkk++; DataOut(kkk)=labITS; kkk++; DataOut(kkk)=lab; kkk++;
1768
1769       for (il=0;il<6;il++) {
1770         idpoint=result.GetIdPoint(il);
1771         idmodule=result.GetIdModule(il);
1772         *(vettid[idmodule]+idpoint)=1;
1773         iotrack->SetIdPoint(il,idpoint);
1774         iotrack->SetIdModule(il,idmodule);
1775       }
1776
1777     //  cout<<"  +++++++++++++  pt e ptg = "<<pt<<" "<<ptg<<"  ++++++++++\n";
1778       Double_t difpt= (pt-ptg)/ptg*100.;
1779       DataOut(kkk)=difpt; kkk++;
1780       Double_t lambdag=TMath::ATan(pzg/ptg);
1781       Double_t   lam=TMath::ATan(tgl);
1782       Double_t diflam = (lam - lambdag)*1000.;
1783       DataOut(kkk) = diflam; kkk++;
1784       Double_t phig=TMath::ATan2(pyg,pxg);  if(phig<0) phig=2.*TMath::Pi()+phig;
1785       Double_t phi=phivertex;
1786
1787       Double_t difphi = (phi - phig)*1000.;
1788       DataOut(kkk)=difphi; kkk++;
1789       DataOut(kkk)=Dtot*1.e4; kkk++;
1790       DataOut(kkk)=Dr*1.e4; kkk++;
1791       DataOut(kkk)=Dz*1.e4; kkk++;
1792       Int_t r;
1793       for (r=0; r<9; r++) { out<<DataOut(r)<<" ";}
1794       out<<"\n";
1795       kkk=0;
1796
1797
1798     } // end if on NumofCluster
1799   //gObjectTable->Print();    // stampa memoria
1800   }  //  end for (int j=min_t; j<=max_t; j++)
1801
1802   out.close();
1803   //outprova.close();
1804
1805   static Bool_t first=kTRUE;
1806   static TFile *tfile;
1807
1808         if(first) {
1809             tfile=new TFile("itstracks.root","RECREATE");
1810             //cout<<"I have opened itstracks.root file "<<endl;
1811         }
1812         first=kFALSE;
1813         tfile->cd();
1814         tfile->ls();
1815
1816    char hname[30];
1817    sprintf(hname,"TreeT%d",evNumber);
1818
1819   tracktree1.Write(hname);
1820
1821
1822
1823             TTree *fAli=gAlice->TreeK();
1824             TFile *fileAli=0;
1825
1826             if (fAli) fileAli =fAli->GetCurrentFile();
1827             fileAli->cd();
1828
1829   ////////////////////////////////////////////////////////////////////////////////////////////////
1830
1831   printf("delete vectors\n");
1832   if(np) delete [] np;
1833   if(vettid) delete [] vettid;
1834
1835 }
1836