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