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