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