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