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