]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCFast.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added pictures and fake calibration
[u/mrichter/AliRoot.git] / TPC / AliTPCFast.cxx
CommitLineData
6d75e4b6 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/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// fast TPC cluster simulation //
21// //
22///////////////////////////////////////////////////////////////////////////////
23
24#include <TParticle.h>
e8d02863 25#include <TVector.h>
6d75e4b6 26
27#include "AliRunLoader.h"
28#include "AliRun.h"
29#include "AliMC.h"
30#include "AliTPC.h"
31#include "AliTPCParam.h"
32#include "AliTPCClustersArray.h"
33#include "AliTPCClustersRow.h"
34#include "AliTPCcluster.h"
35#include "AliComplexCluster.h"
36#include "AliTPCFast.h"
8c2b3fd7 37#include "AliLog.h"
6d75e4b6 38
39ClassImp(AliTPCFast)
40
41
42//_____________________________________________________________________________
43void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
44{
45 //--------------------------------------------------------
46 // TPC simple cluster generator from hits
47 // obtained from the TPC Fast Simulator
48 // The point errors are taken from the parametrization
49 //--------------------------------------------------------
50
51 //-----------------------------------------------------------------
52 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
53 //-----------------------------------------------------------------
54 // Adopted to Marian's cluster data structure by I.Belikov, CERN,
55 // Jouri.Belikov@cern.ch
56 //----------------------------------------------------------------
57
58 /////////////////////////////////////////////////////////////////////////////
59 //
60 //---------------------------------------------------------------------
61 // ALICE TPC Cluster Parameters
62 //--------------------------------------------------------------------
63
64
65
66 // Cluster width in rphi
67 const Float_t kACrphi=0.18322;
68 const Float_t kBCrphi=0.59551e-3;
69 const Float_t kCCrphi=0.60952e-1;
70 // Cluster width in z
71 const Float_t kACz=0.19081;
72 const Float_t kBCz=0.55938e-3;
73 const Float_t kCCz=0.30428;
74
75
76 AliLoader* loader = runLoader->GetLoader("TPCLoader");
77 if (!loader) {
8c2b3fd7 78 AliError("No TPC loader found");
6d75e4b6 79 return;
80 }
81 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
82 AliRun* aliRun = runLoader->GetAliRun();
83 if (!aliRun) {
8c2b3fd7 84 AliError("Couldn't get AliRun object");
6d75e4b6 85 return;
86 }
87 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
88 if (!tpc) {
8c2b3fd7 89 AliError("Couldn't get TPC detector");
6d75e4b6 90 return;
91 }
92 AliTPCParam* param = tpc->GetParam();
93 if (!param) {
8c2b3fd7 94 AliError("No TPC parameters available");
6d75e4b6 95 return;
96 }
97
98 //if(fDefaults == 0) SetDefaults();
99
100 Float_t sigmaRphi,sigmaZ,clRphi,clZ;
101 //
102 TParticle *particle; // pointer to a given particle
103 AliTPChit *tpcHit; // pointer to a sigle TPC hit
104 Int_t sector;
105 Int_t ipart;
106 Float_t xyz[5];
107 Float_t pl,pt,tanth,rpad,ratio;
108 Float_t cph,sph;
109
110 //---------------------------------------------------------------
111 // Get the access to the tracks
112 //---------------------------------------------------------------
113
114 TTree *tH = loader->TreeH();
115 if (tH == 0x0)
8c2b3fd7 116 AliFatal("Can not find TreeH in folder");
117
6d75e4b6 118 tpc->SetTreeAddress();
119
120 Stat_t ntracks = tH->GetEntries();
121
122 //Switch to the output file
123
124 if (loader->TreeR() == 0x0) loader->MakeTree("R");
125
8c2b3fd7 126 AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
6d75e4b6 127
128 runLoader->CdGAFile();
129 //param->Write(param->GetTitle());
130
131 AliTPCClustersArray carray;
132 carray.Setup(param);
133 carray.SetClusterType("AliTPCcluster");
134 carray.MakeTree(loader->TreeR());
135
136 Int_t nclusters=0; //cluster counter
137
138 //------------------------------------------------------------
139 // Loop over all sectors (72 sectors for 20 deg
140 // segmentation for both lower and upper sectors)
141 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
142 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
143 //
144 // First cluster for sector 0 starts at "0"
145 //------------------------------------------------------------
146
147 for(Int_t isec=0;isec<param->GetNSector();isec++){
148 //MI change
149 param->AdjustCosSin(isec,cph,sph);
150
151 //------------------------------------------------------------
152 // Loop over tracks
153 //------------------------------------------------------------
154
155 for(Int_t track=0;track<ntracks;track++){
156 tpc->ResetHits();
157 tH->GetEvent(track);
158 //
159 // Get number of the TPC hits
160 //
161 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
162
163 // Loop over hits
164 //
165 while(tpcHit){
8ebf1d48 166
6d75e4b6 167 if (tpcHit->fQ == 0.) {
168 tpcHit = (AliTPChit*) tpc->NextHit();
169 continue; //information about track (I.Belikov)
170 }
8ebf1d48 171 sector=tpcHit->fSector; // sector number
172
173 if(sector != isec){
174 tpcHit = (AliTPChit*) tpc->NextHit();
175 continue;
176 }
177 ipart=tpcHit->Track();
178 particle=aliRun->GetMCApp()->Particle(ipart);
179 pl=particle->Pz();
180 pt=particle->Pt();
181 if(pt < 1.e-9) pt=1.e-9;
182 tanth=pl/pt;
183 tanth = TMath::Abs(tanth);
184 rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
185 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
186
187 // space-point resolutions
188
6d75e4b6 189 sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
190 sigmaZ =AliTPCcluster::SigmaZ2(rpad,tanth );
191
192 // cluster widths
193
194 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
195 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
196
197 // temporary protection
198
199 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
200 if(sigmaZ < 0.) sigmaZ=0.4e-3;
201 if(clRphi < 0.) clRphi=2.5e-3;
202 if(clZ < 0.) clZ=2.5e-5;
203
204 //
205
206 //
207 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
208 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
209 //
210 Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
211 Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
212 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
213 Float_t alpha=(isec < param->GetNInnerSector()) ?
214 param->GetInnerAngle() : param->GetOuterAngle();
215 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
216 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
217 xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
218 if (TMath::Abs(xyz[1])>param->GetZLength()) xyz[1]=tpcHit->Z();
219 xyz[2]=sigmaRphi; // fSigmaY2
220 xyz[3]=sigmaZ; // fSigmaZ2
221 xyz[4]=tpcHit->fQ; // q
222
223 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
224 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
225
226 Int_t tracks[3]={tpcHit->Track(), -1, -1};
227 AliTPCcluster cluster(tracks,xyz);
228
229 clrow->InsertCluster(&cluster); nclusters++;
230
231 tpcHit = (AliTPChit*)tpc->NextHit();
232
233
234 } // end of loop over hits
235
236 } // end of loop over tracks
237
238 Int_t nrows=param->GetNRow(isec);
239 for (Int_t irow=0; irow<nrows; irow++) {
240 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
241 if (!clrow) continue;
242 carray.StoreRow(isec,irow);
243 carray.ClearRow(isec,irow);
244 }
245
246 } // end of loop over sectors
247
248 // cerr<<"Number of made clusters : "<<nclusters<<" \n";
249 loader->WriteRecPoints("OVERWRITE");
250
251
252} // end of function
253
8ebf1d48 254
255
256//_________________________________________________________________
257
258
259
260void AliTPCFast::Hits2ExactClusters(AliRunLoader* runLoader) const{
261
262
263
264 AliLoader* loader = runLoader->GetLoader("TPCLoader");
265 if (!loader) {
266 AliError("No TPC loader found");
267 return;
268 }
269 AliTPC* tpc = (AliTPC*) gAlice->GetDetector("TPC");
270 if (!tpc) {
271 AliError("Couldn't get TPC detector");
272 return;
273 }
274 AliTPCParam* param = tpc->GetParam();
275 if (!param) {
276 AliError("No TPC parameters available");
277 return;
278 }
279 //---------------------------------------------------------------
280 // Get the access to the tracks
281 //---------------------------------------------------------------
282
283 TTree *tH = loader->TreeH();
284 if (tH == 0x0) { runLoader->LoadHits("TPC"); loader->TreeH();}
285
286 tpc->SetTreeAddress();
287
288
289 //Switch to the output file
290
291 if (loader->TreeR() == 0x0) loader->MakeTree("R");
292
293 AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
294
295 runLoader->CdGAFile();
296 //param->Write(param->GetTitle());
297
298 AliTPCClustersArray carray;
299 carray.Setup(param);
300 carray.SetClusterType("AliTPCclusterMI");
301 carray.MakeTree(loader->TreeR());
302
303 //------------------------------------------------------------
304 // Loop over all sectors (72 sectors for 20 deg)
305 //------------------------------------------------------------
306
307 for(Int_t isec=0;isec<param->GetNSector();isec++){
308 Hits2ExactClustersSector(runLoader, &carray, isec);
309 } // end of loop over sectors
310 loader->WriteRecPoints("OVERWRITE");
311}
312
313
314
315
6d75e4b6 316//_________________________________________________________________
317void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
318 AliTPCClustersArray* clustersArray,
319 Int_t isec) const
320{
321 //--------------------------------------------------------
322 //calculate exact cross point of track and given pad row
8ebf1d48 323 //resulting values are expressed in "local" coordinata
324 //the sigmaY2 and sigmaZ2 of cluster are the shape of cluster parameters
325 // - thay are used later on for error parameterization
6d75e4b6 326 //--------------------------------------------------------
327
328 //-----------------------------------------------------------------
329 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
330 //-----------------------------------------------------------------
331 //
332 if (clustersArray==0){
333 return;
334 }
335 AliLoader* loader = runLoader->GetLoader("TPCLoader");
336 if (!loader) {
8c2b3fd7 337 AliError("No TPC loader found");
6d75e4b6 338 return;
339 }
340 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
341 AliRun* aliRun = runLoader->GetAliRun();
342 if (!aliRun) {
8c2b3fd7 343 AliError("Couldn't get AliRun object");
6d75e4b6 344 return;
345 }
346 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
347 if (!tpc) {
8c2b3fd7 348 AliError("Couldn't get TPC detector");
6d75e4b6 349 return;
350 }
351 AliTPCParam* param = tpc->GetParam();
352 if (!param) {
8c2b3fd7 353 AliError("No TPC parameters available");
6d75e4b6 354 return;
355 }
356 //
8ebf1d48 357 //
6d75e4b6 358 AliTPChit *tpcHit; // pointer to a sigle TPC hit
359 // Int_t sector,nhits;
360 Int_t ipart;
361 const Int_t kcmaxhits=30000;
8c2b3fd7 362 TVector * xxxx = new TVector(kcmaxhits*4);
363 TVector & xxx = *xxxx;
6d75e4b6 364 Int_t maxhits = kcmaxhits;
365 //construct array for each padrow
366 for (Int_t i=0; i<param->GetNRow(isec);i++)
367 clustersArray->CreateRow(isec,i);
368
369 //---------------------------------------------------------------
370 // Get the access to the tracks
371 //---------------------------------------------------------------
372
373 TTree *tH = loader->TreeH();
374 if (tH == 0x0)
8c2b3fd7 375 AliFatal("Can not find TreeH in folder");
376
6d75e4b6 377 tpc->SetTreeAddress();
378
379 Stat_t ntracks = tH->GetEntries();
6d75e4b6 380 //MI change
381 TBranch * branch=0;
382 if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
383 else branch = tH->GetBranch("TPC");
384
385 //------------------------------------------------------------
386 // Loop over tracks
387 //------------------------------------------------------------
388
389 for(Int_t track=0;track<ntracks;track++){
390 Bool_t isInSector=kTRUE;
391 tpc->ResetHits();
8ebf1d48 392 isInSector = tpc->TrackInVolume(isec,track);
6d75e4b6 393 if (!isInSector) continue;
394 //MI change
395 branch->GetEntry(track); // get next track
396 //
6d75e4b6 397 // Loop over hits
398 //
399 Int_t currentIndex=0;
400 Int_t lastrow=-1; //last writen row
401
402 //M.I. changes
403
404 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
405 while(tpcHit){
406
407 Int_t sector=tpcHit->fSector; // sector number
408 if(sector != isec){
409 tpcHit = (AliTPChit*) tpc->NextHit();
410 continue;
411 }
412
413 ipart=tpcHit->Track();
6d75e4b6 414
415 //find row number
416
417 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
418 Int_t index[3]={1,isec,0};
8ebf1d48 419 Int_t currentrow = param->GetPadRow(x,index);
420
6d75e4b6 421 if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
422 if (lastrow<0) lastrow=currentrow;
8ebf1d48 423 //
424 //
425 if (currentrow!=lastrow){
6d75e4b6 426 if (currentIndex>2){
427 Float_t sumx=0;
428 Float_t sumx2=0;
429 Float_t sumx3=0;
430 Float_t sumx4=0;
431 Float_t sumy=0;
432 Float_t sumxy=0;
433 Float_t sumx2y=0;
434 Float_t sumz=0;
435 Float_t sumxz=0;
436 Float_t sumx2z=0;
437 Float_t sumq=0;
438 for (Int_t index=0;index<currentIndex;index++){
439 Float_t x,x2,x3,x4;
440 x=x2=x3=x4=xxx(index*4);
441 x2*=x;
442 x3*=x2;
443 x4*=x3;
444 sumx+=x;
445 sumx2+=x2;
446 sumx3+=x3;
447 sumx4+=x4;
448 sumy+=xxx(index*4+1);
449 sumxy+=xxx(index*4+1)*x;
450 sumx2y+=xxx(index*4+1)*x2;
451 sumz+=xxx(index*4+2);
452 sumxz+=xxx(index*4+2)*x;
453 sumx2z+=xxx(index*4+2)*x2;
454 sumq+=xxx(index*4+3);
455 }
6d75e4b6 456 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
457 sumx2*(sumx*sumx3-sumx2*sumx2);
458
459 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
460 sumx2*(sumxy*sumx3-sumx2y*sumx2);
461 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
462 sumx2*(sumxz*sumx3-sumx2z*sumx2);
463
464 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
465 sumx2*(sumx*sumx2y-sumx2*sumxy);
466 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
467 sumx2*(sumx*sumx2z-sumx2*sumxz);
468
8ebf1d48 469 if (TMath::Abs(det)<0.00000000000000001){
6d75e4b6 470 tpcHit = (AliTPChit*)tpc->NextHit();
471 continue;
472 }
473
8ebf1d48 474 Float_t y=detay/det;
475 Float_t z=detaz/det;
6d75e4b6 476 Float_t by=detby/det; //y angle
477 Float_t bz=detbz/det; //z angle
478 sumy/=Float_t(currentIndex);
479 sumz/=Float_t(currentIndex);
480
6d75e4b6 481
8ebf1d48 482 //
483 //
484 Float_t sign = (tpcHit->Z()>0)? 1.:-1.;
485 y = (y+0.5)*param->GetPadPitchWidth(isec);
486 z = z*param->GetZWidth();
487 //
488 // y expected shape
489 Double_t sigmay2 = z*fParam->GetDiffL();
490 sigmay2 *= sigmay2;
491 Float_t angulary = by*param->GetPadPitchLength(isec,lastrow);
492 angulary*=angulary;
493 angulary/=12;
494 sigmay2 +=angulary+0.25*param->GetPadPitchWidth(isec)*param->GetPadPitchWidth(isec);
495 //
496 // z expected shape
497 Double_t sigmaz2 = z*fParam->GetDiffT();
498 sigmaz2 *= sigmaz2;
499 Float_t angularz = bz*param->GetPadPitchLength(isec,lastrow);
500 angularz*=angularz;
501 angularz/=12;
502 sigmaz2 +=angularz+0.25;
503 //
504 sigmaz2 = TMath::Min(sigmaz2,1.);
505 sigmay2 = TMath::Min(sigmay2,1.);
506 //
507 //
508 z = sign*(param->GetZLength() - z);
509 if (TMath::Abs(z)< param->GetZLength()-1){
510 AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
511 if (row!=0) {
512 AliTPCclusterMI* cl = new((AliTPCclusterMI*)row->Append()) AliTPCclusterMI ;
513 // AliTPCclusterMI cl;
514 cl->SetZ(z);
515 cl->SetY(y);
516 cl->SetQ(sumq);
517 if (TMath::Abs(sumq)<1000){
518 cl->SetMax(Short_t(sumq));
519 }
520 else{
521 cl->SetMax(0);
522 }
523 cl->SetSigmaZ2(sigmaz2);
524 cl->SetSigmaY2(sigmay2);
525 cl->SetLabel(ipart,0);
526 cl->SetLabel(-1,1);
527 cl->SetLabel(-1,2);
528 cl->SetType(0);
529 }
530 }
531 } //end of calculating cluster for given row
532 currentIndex=0;
533 lastrow=currentrow;
534 } // end of crossing rows
535 //
6d75e4b6 536 if ( currentIndex>=maxhits){
537 maxhits+=kcmaxhits;
538 xxx.ResizeTo(4*maxhits);
539 }
540 xxx(currentIndex*4)=x[0];
541 xxx(currentIndex*4+1)=x[1];
542 xxx(currentIndex*4+2)=x[2];
543 xxx(currentIndex*4+3)=tpcHit->fQ;
544 currentIndex++;
8ebf1d48 545 tpcHit = (AliTPChit*)tpc->NextHit();
546 } // end of loop over hits
547 } // end of loop over tracks
548 //write padrows to tree
549 for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
550 clustersArray->StoreRow(isec,ii);
551 clustersArray->ClearRow(isec,ii);
552 }
6d75e4b6 553 xxxx->Delete();
8ebf1d48 554
6d75e4b6 555}
556
8ebf1d48 557