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