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