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