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