1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // fast TPC cluster simulation //
22 ///////////////////////////////////////////////////////////////////////////////
24 #include <TParticle.h>
26 #include "AliRunLoader.h"
30 #include "AliTPCParam.h"
31 #include "AliTPCClustersArray.h"
32 #include "AliTPCClustersRow.h"
33 #include "AliTPCcluster.h"
34 #include "AliComplexCluster.h"
35 #include "AliTPCFast.h"
41 //_____________________________________________________________________________
42 void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
44 //--------------------------------------------------------
45 // TPC simple cluster generator from hits
46 // obtained from the TPC Fast Simulator
47 // The point errors are taken from the parametrization
48 //--------------------------------------------------------
50 //-----------------------------------------------------------------
51 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
52 //-----------------------------------------------------------------
53 // Adopted to Marian's cluster data structure by I.Belikov, CERN,
54 // Jouri.Belikov@cern.ch
55 //----------------------------------------------------------------
57 /////////////////////////////////////////////////////////////////////////////
59 //---------------------------------------------------------------------
60 // ALICE TPC Cluster Parameters
61 //--------------------------------------------------------------------
65 // Cluster width in rphi
66 const Float_t kACrphi=0.18322;
67 const Float_t kBCrphi=0.59551e-3;
68 const Float_t kCCrphi=0.60952e-1;
70 const Float_t kACz=0.19081;
71 const Float_t kBCz=0.55938e-3;
72 const Float_t kCCz=0.30428;
75 AliLoader* loader = runLoader->GetLoader("TPCLoader");
77 Error("Hits2Clusters", "No TPC loader found");
80 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
81 AliRun* aliRun = runLoader->GetAliRun();
83 Error("Hits2Clusters", "Couldn't get AliRun object");
86 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
88 Error("Hits2Clusters", "Couldn't get TPC detector");
91 AliTPCParam* param = tpc->GetParam();
93 Error("Hits2Clusters", "No TPC parameters available");
97 //if(fDefaults == 0) SetDefaults();
99 Float_t sigmaRphi,sigmaZ,clRphi,clZ;
101 TParticle *particle; // pointer to a given particle
102 AliTPChit *tpcHit; // pointer to a sigle TPC hit
106 Float_t pl,pt,tanth,rpad,ratio;
109 //---------------------------------------------------------------
110 // Get the access to the tracks
111 //---------------------------------------------------------------
113 TTree *tH = loader->TreeH();
116 Fatal("Hits2Clusters","Can not find TreeH in folder");
119 tpc->SetTreeAddress();
121 Stat_t ntracks = tH->GetEntries();
123 //Switch to the output file
125 if (loader->TreeR() == 0x0) loader->MakeTree("R");
127 cout<<"param->GetTitle() = "<<param->GetTitle()<<endl;
129 runLoader->CdGAFile();
130 //param->Write(param->GetTitle());
132 AliTPCClustersArray carray;
134 carray.SetClusterType("AliTPCcluster");
135 carray.MakeTree(loader->TreeR());
137 Int_t nclusters=0; //cluster counter
139 //------------------------------------------------------------
140 // Loop over all sectors (72 sectors for 20 deg
141 // segmentation for both lower and upper sectors)
142 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
143 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
145 // First cluster for sector 0 starts at "0"
146 //------------------------------------------------------------
148 for(Int_t isec=0;isec<param->GetNSector();isec++){
150 param->AdjustCosSin(isec,cph,sph);
152 //------------------------------------------------------------
154 //------------------------------------------------------------
156 for(Int_t track=0;track<ntracks;track++){
160 // Get number of the TPC hits
162 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
168 if (tpcHit->fQ == 0.) {
169 tpcHit = (AliTPChit*) tpc->NextHit();
170 continue; //information about track (I.Belikov)
172 sector=tpcHit->fSector; // sector number
175 tpcHit = (AliTPChit*) tpc->NextHit();
178 ipart=tpcHit->Track();
179 particle=aliRun->GetMCApp()->Particle(ipart);
182 if(pt < 1.e-9) pt=1.e-9;
184 tanth = TMath::Abs(tanth);
185 rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
186 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
188 // space-point resolutions
190 sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
191 sigmaZ =AliTPCcluster::SigmaZ2(rpad,tanth );
195 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
196 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
198 // temporary protection
200 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
201 if(sigmaZ < 0.) sigmaZ=0.4e-3;
202 if(clRphi < 0.) clRphi=2.5e-3;
203 if(clZ < 0.) clZ=2.5e-5;
208 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
209 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
211 Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
212 Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
213 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
214 Float_t alpha=(isec < param->GetNInnerSector()) ?
215 param->GetInnerAngle() : param->GetOuterAngle();
216 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
217 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
218 xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
219 if (TMath::Abs(xyz[1])>param->GetZLength()) xyz[1]=tpcHit->Z();
220 xyz[2]=sigmaRphi; // fSigmaY2
221 xyz[3]=sigmaZ; // fSigmaZ2
222 xyz[4]=tpcHit->fQ; // q
224 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
225 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
227 Int_t tracks[3]={tpcHit->Track(), -1, -1};
228 AliTPCcluster cluster(tracks,xyz);
230 clrow->InsertCluster(&cluster); nclusters++;
232 tpcHit = (AliTPChit*)tpc->NextHit();
235 } // end of loop over hits
237 } // end of loop over tracks
239 Int_t nrows=param->GetNRow(isec);
240 for (Int_t irow=0; irow<nrows; irow++) {
241 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
242 if (!clrow) continue;
243 carray.StoreRow(isec,irow);
244 carray.ClearRow(isec,irow);
247 } // end of loop over sectors
249 // cerr<<"Number of made clusters : "<<nclusters<<" \n";
250 loader->WriteRecPoints("OVERWRITE");
255 //_________________________________________________________________
256 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
257 AliTPCClustersArray* clustersArray,
260 //--------------------------------------------------------
261 //calculate exact cross point of track and given pad row
262 //resulting values are expressed in "digit" coordinata
263 //--------------------------------------------------------
265 //-----------------------------------------------------------------
266 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
267 //-----------------------------------------------------------------
269 if (clustersArray==0){
272 AliLoader* loader = runLoader->GetLoader("TPCLoader");
274 Error("Hits2ExactClustersSector", "No TPC loader found");
277 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
278 AliRun* aliRun = runLoader->GetAliRun();
280 Error("Hits2ExactClustersSector", "Couldn't get AliRun object");
283 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
285 Error("Hits2ExactClustersSector", "Couldn't get TPC detector");
288 AliTPCParam* param = tpc->GetParam();
290 Error("Hits2ExactClustersSector", "No TPC parameters available");
294 TParticle *particle; // pointer to a given particle
295 AliTPChit *tpcHit; // pointer to a sigle TPC hit
296 // Int_t sector,nhits;
298 const Int_t kcmaxhits=30000;
299 AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
300 AliTPCFastVector & xxx = *xxxx;
301 Int_t maxhits = kcmaxhits;
302 //construct array for each padrow
303 for (Int_t i=0; i<param->GetNRow(isec);i++)
304 clustersArray->CreateRow(isec,i);
306 //---------------------------------------------------------------
307 // Get the access to the tracks
308 //---------------------------------------------------------------
310 TTree *tH = loader->TreeH();
313 Fatal("Hits2Clusters","Can not find TreeH in folder");
316 tpc->SetTreeAddress();
318 Stat_t ntracks = tH->GetEntries();
319 Int_t npart = aliRun->GetMCApp()->GetNtrack();
322 if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
323 else branch = tH->GetBranch("TPC");
325 //------------------------------------------------------------
327 //------------------------------------------------------------
329 for(Int_t track=0;track<ntracks;track++){
330 Bool_t isInSector=kTRUE;
332 isInSector = tpc->TrackInVolume(isec,track);
333 if (!isInSector) continue;
335 branch->GetEntry(track); // get next track
337 // Get number of the TPC hits and a pointer
341 Int_t currentIndex=0;
342 Int_t lastrow=-1; //last writen row
346 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
349 Int_t sector=tpcHit->fSector; // sector number
351 tpcHit = (AliTPChit*) tpc->NextHit();
355 ipart=tpcHit->Track();
356 if (ipart<npart) particle=aliRun->GetMCApp()->Particle(ipart);
360 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
361 Int_t index[3]={1,isec,0};
362 Int_t currentrow = param->GetPadRow(x,index) ;
363 if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
364 if (lastrow<0) lastrow=currentrow;
365 if (currentrow==lastrow){
366 if ( currentIndex>=maxhits){
368 xxx.ResizeTo(4*maxhits);
370 xxx(currentIndex*4)=x[0];
371 xxx(currentIndex*4+1)=x[1];
372 xxx(currentIndex*4+2)=x[2];
373 xxx(currentIndex*4+3)=tpcHit->fQ;
389 for (Int_t index=0;index<currentIndex;index++){
391 x=x2=x3=x4=xxx(index*4);
399 sumy+=xxx(index*4+1);
400 sumxy+=xxx(index*4+1)*x;
401 sumx2y+=xxx(index*4+1)*x2;
402 sumz+=xxx(index*4+2);
403 sumxz+=xxx(index*4+2)*x;
404 sumx2z+=xxx(index*4+2)*x2;
405 sumq+=xxx(index*4+3);
407 Float_t centralPad = (param->GetNPads(isec,lastrow)-1)/2;
408 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
409 sumx2*(sumx*sumx3-sumx2*sumx2);
411 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
412 sumx2*(sumxy*sumx3-sumx2y*sumx2);
413 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
414 sumx2*(sumxz*sumx3-sumx2z*sumx2);
416 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
417 sumx2*(sumx*sumx2y-sumx2*sumxy);
418 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
419 sumx2*(sumx*sumx2z-sumx2*sumxz);
421 if (TMath::Abs(det)<0.00001){
422 tpcHit = (AliTPChit*)tpc->NextHit();
426 Float_t y=detay/det+centralPad;
428 Float_t by=detby/det; //y angle
429 Float_t bz=detbz/det; //z angle
430 sumy/=Float_t(currentIndex);
431 sumz/=Float_t(currentIndex);
433 AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
435 AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
436 // AliComplexCluster cl;
442 cl->fTracks[0]=ipart;
446 } //end of calculating cluster for given row
449 tpcHit = (AliTPChit*)tpc->NextHit();
450 } // end of loop over hits
451 } // end of loop over tracks
452 //write padrows to tree
453 for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
454 clustersArray->StoreRow(isec,ii);
455 clustersArray->ClearRow(isec,ii);
462 void AliTPCFast::FindTrackHitsIntersection(AliRunLoader* runLoader,
463 AliTPCClustersArray* clustersArray) const
467 //fill clones array with intersection of current point with the
469 AliLoader* loader = runLoader->GetLoader("TPCLoader");
471 Error("FindTrackHitsIntersection", "No TPC loader found");
474 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
475 AliRun* aliRun = runLoader->GetAliRun();
477 Error("FindTrackHitsIntersection", "Couldn't get AliRun object");
480 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
482 Error("FindTrackHitsIntersection", "Couldn't get TPC detector");
485 AliTPCParam* param = tpc->GetParam();
487 Error("FindTrackHitsIntersection", "No TPC parameters available");
494 const Int_t kcmaxhits=30000;
495 AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
496 AliTPCFastVector & xxx = *xxxx;
497 Int_t maxhits = kcmaxhits;
500 AliTPChit * tpcHit=0;
501 tpcHit = (AliTPChit*)tpc->FirstHit2(-1);
502 Int_t currentIndex=0;
503 Int_t lastrow=-1; //last writen row
506 if (tpcHit==0) continue;
507 sector=tpcHit->fSector; // sector number
508 ipart=tpcHit->Track();
512 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
513 Int_t index[3]={1,sector,0};
514 Int_t currentrow = param->GetPadRow(x,index) ;
515 if (currentrow<0) continue;
516 if (lastrow<0) lastrow=currentrow;
517 if (currentrow==lastrow){
518 if ( currentIndex>=maxhits){
520 xxx.ResizeTo(4*maxhits);
522 xxx(currentIndex*4)=x[0];
523 xxx(currentIndex*4+1)=x[1];
524 xxx(currentIndex*4+2)=x[2];
525 xxx(currentIndex*4+3)=tpcHit->fQ;
541 for (Int_t index=0;index<currentIndex;index++){
543 x=x2=x3=x4=xxx(index*4);
551 sumy+=xxx(index*4+1);
552 sumxy+=xxx(index*4+1)*x;
553 sumx2y+=xxx(index*4+1)*x2;
554 sumz+=xxx(index*4+2);
555 sumxz+=xxx(index*4+2)*x;
556 sumx2z+=xxx(index*4+2)*x2;
557 sumq+=xxx(index*4+3);
559 Float_t centralPad = (param->GetNPads(sector,lastrow)-1)/2;
560 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
561 sumx2*(sumx*sumx3-sumx2*sumx2);
563 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
564 sumx2*(sumxy*sumx3-sumx2y*sumx2);
565 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
566 sumx2*(sumxz*sumx3-sumx2z*sumx2);
568 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
569 sumx2*(sumx*sumx2y-sumx2*sumxy);
570 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
571 sumx2*(sumx*sumx2z-sumx2*sumxz);
573 Float_t y=detay/det+centralPad;
575 Float_t by=detby/det; //y angle
576 Float_t bz=detbz/det; //z angle
577 sumy/=Float_t(currentIndex);
578 sumz/=Float_t(currentIndex);
580 AliComplexCluster cl;
588 AliTPCClustersRow * row = (clustersArray->GetRow(sector,lastrow));
589 if (row!=0) row->InsertCluster(&cl);
592 } //end of calculating cluster for given row
594 } // end of loop over hits