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