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