Transition to NewIO
[u/mrichter/AliRoot.git] / ITS / AliITSFindClustersV2.C
CommitLineData
88cb7938 1/****************************************************************
2* This macro converts AliITSRecPoint(s) to AliITSclusterV2(s) *
3* Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch *
4*****************************************************************/
f2718ade 5
03248e6d 6#ifndef __CINT__
88cb7938 7 #include <Riostream.h>
03248e6d 8
9 #include "AliRun.h"
88cb7938 10 #include "AliRunLoader.h"
11 #include "AliLoader.h"
03248e6d 12 #include "AliITS.h"
13 #include "AliITSgeom.h"
f2718ade 14 #include "AliITSclustererV2.h"
03248e6d 15
f2718ade 16 #include "TStopwatch.h"
03248e6d 17#endif
18
88cb7938 19Int_t AliITSFindClustersV2(Char_t SlowOrFast='f')
20{
03248e6d 21
88cb7938 22 cerr<<"AliITSRecPoint(s) -> AliITSclusterV2(s)...\n";
23
24 if (gAlice)
25 {
26 delete gAlice->GetRunLoader();
27 delete gAlice;
28 gAlice=0;
29 }
30
31 AliRunLoader* rl = AliRunLoader::Open("galice.root");
32 if (rl == 0x0)
33 {
34 cerr<<"AliITSHits2DigitsDefault.C : Can not open session RL=NULL"
35 << endl;
36 return 3;
37 }
38
39 Int_t retval = rl->LoadgAlice();
40 if (retval)
41 {
42 cerr<<"AliITSHits2DigitsDefault.C : LoadgAlice returned error"
43 << endl;
44 delete rl;
45 return 3;
46 }
47 gAlice=rl->GetAliRun();
48 rl->LoadHeader();
49 retval = rl->LoadKinematics();
50 if (retval)
51 {
52 cerr<<"AliITSHits2DigitsDefault.C : LoadKinematics returned error"
53 << endl;
54 delete rl;
55 return 3;
56 }
57
58 AliITSLoader* gime = (AliITSLoader*)rl->GetLoader("ITSLoader");
59 if (gime == 0x0)
60 {
61 cerr<<"AliITSHits2DigitsDefault.C : can not get ITS loader"
62 << endl;
63 }
03248e6d 64
88cb7938 65 rl->GetEvent(0);
f2718ade 66
88cb7938 67 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
68 if (!ITS) { cerr<<"Can't find the ITS !\n"; delete rl; return 3; }
03248e6d 69 AliITSgeom *geom=ITS->GetITSgeom();
88cb7938 70
71 TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
72
73 gime->LoadRawClusters("recreate");
f2718ade 74
88cb7938 75 if (SlowOrFast=='f')
76 {
77 gime->SetRecPointsFileName("ITS.FastRecPoints.root");
78 }
79 if (gime->LoadRecPoints())
80 {
81 cerr<<"Load Rec Pints returned error !\n";
82 delete rl;
83 return 4;
84 }
85
86 TClonesArray *points = new TClonesArray("AliITSRecPoint",10000);
87
88 Float_t lp[5];
89 Int_t lab[6];
90
91 Int_t iEvent;
92 for (iEvent = 0; iEvent< rl->GetNumberOfEvents() ; iEvent++)
93 {
94
95 rl->GetEvent(iEvent);
96
97 TTree *cTree = gime->TreeC();
98 if (cTree == 0x0)
99 {
100 gime->MakeTree("C");
101 cTree = gime->TreeC();
102 }
103
104 cTree->Branch("Clusters",&clusters);
105
106 TTree *pTree=gime->TreeR();
107 if (pTree == 0x0)
108 {
109 cerr<<"Can not get TreeR !\n"; delete rl;return 5;
110 }
111 TBranch *branch = 0;
112 if (SlowOrFast=='f') {
113 branch = pTree->GetBranch("ITSRecPointsF");
114 }
115 else {
116 branch = pTree->GetBranch("ITSRecPoints");
117 }
118 if (!branch)
119 {
120 cerr<<"Can't get ITSRecPoints branch !\n";
121 delete rl;
122 return 6;
123 }
124
125 branch->SetAddress(&points);
126
127 AliStack* stack = rl->Stack();
128 if (stack == 0x0)
129 {
130 cerr<<"AliITSFindClustersV2.C : Can not get stack"
131 << endl;
132 delete rl;
133 return 3;
134 }
135
136 TClonesArray &cl=*clusters;
137 Int_t nclusters=0;
138 Int_t nentr=(Int_t)branch->GetEntries();
139
140 cerr<<"Number of entries: "<<nentr<<endl;
141
142 for (Int_t i=0; i<nentr; i++)
143 {
144 points->Clear();
145 clusters->Clear();
146 branch->GetEvent(i);
147 Int_t ncl=points->GetEntriesFast(); if (ncl==0){cTree->Fill();continue;}
148 Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
149 if ( (lay<0) || (lad<0) || (det<0))
150 {
151 ::Error("AliITSFindClustersV2.C","No such a module %d",i);
152 continue;
153 }
154 Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
155 Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot);
156 Double_t yshift = x*rot[0] + y*rot[1];
157 Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
158 nclusters+=ncl;
159
160 Float_t kmip=1; // ADC->mip normalization factor for the SDD and SSD
161 if(lay==4 || lay==3){kmip=280.;};
162 if(lay==6 || lay==5){kmip=38.;};
163
164 for (Int_t j=0; j<ncl; j++)
165 {
166 AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
167 //Float_t lp[5];
168 lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
169 lp[1]=p->GetZ()+zshift;
170 lp[2]=p->GetSigmaX2();
171 lp[3]=p->GetSigmaZ2();
172 lp[4]=p->GetQ(); lp[4]/=kmip;
173 //Int_t lab[6];
174 lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
175 lab[3]=ndet;
03248e6d 176
88cb7938 177 Int_t label=lab[0];
178 if (label>=0) {
179 TParticle *part=(TParticle*)stack->Particle(label);
180 if (part == 0x0)
181 cerr<<"Can not get particle with label "<<label<<endl;
182 label=-3;
183 while (part->P() < 0.005) {
184 Int_t m=part->GetFirstMother();
185 if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
186 label=m;
187 part=(TParticle*)gAlice->Particle(label);
188 }
189 if (lab[1]<0) lab[1]=label;
190 else if (lab[2]<0) lab[2]=label;
191 else cerr<<"No empty labels !\n";
192 }
03248e6d 193
88cb7938 194 new(cl[j]) AliITSclusterV2(lab,lp);
195 }
196 cTree->Fill();
197// clusters->Delete(); points->Delete();
198 }
199 gime->WriteRawClusters("OVERWRITE");
200 cerr<<"Number of clusters: "<<nclusters<<endl;
201
202 }
203
204 delete clusters;
205 delete points;
206 delete rl;
03248e6d 207 return 0;
88cb7938 208
03248e6d 209}
210
211
212
88cb7938 213
214
215
216
217
218
219
220
221
222
223