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 AliError("No TPC loader found");
80 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
81 AliRun* aliRun = runLoader->GetAliRun();
83 AliError("Couldn't get AliRun object");
86 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
88 AliError("Couldn't get TPC detector");
91 AliTPCParam* param = tpc->GetParam();
93 AliError("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();
115 AliFatal("Can not find TreeH in folder");
117 tpc->SetTreeAddress();
119 Stat_t ntracks = tH->GetEntries();
121 //Switch to the output file
123 if (loader->TreeR() == 0x0) loader->MakeTree("R");
125 AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
127 runLoader->CdGAFile();
128 //param->Write(param->GetTitle());
130 AliTPCClustersArray carray;
132 carray.SetClusterType("AliTPCcluster");
133 carray.MakeTree(loader->TreeR());
135 Int_t nclusters=0; //cluster counter
137 //------------------------------------------------------------
138 // Loop over all sectors (72 sectors for 20 deg
139 // segmentation for both lower and upper sectors)
140 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
141 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
143 // First cluster for sector 0 starts at "0"
144 //------------------------------------------------------------
146 for(Int_t isec=0;isec<param->GetNSector();isec++){
148 param->AdjustCosSin(isec,cph,sph);
150 //------------------------------------------------------------
152 //------------------------------------------------------------
154 for(Int_t track=0;track<ntracks;track++){
158 // Get number of the TPC hits
160 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
166 if (tpcHit->fQ == 0.) {
167 tpcHit = (AliTPChit*) tpc->NextHit();
168 continue; //information about track (I.Belikov)
170 sector=tpcHit->fSector; // sector number
173 tpcHit = (AliTPChit*) tpc->NextHit();
176 ipart=tpcHit->Track();
177 particle=aliRun->GetMCApp()->Particle(ipart);
180 if(pt < 1.e-9) pt=1.e-9;
182 tanth = TMath::Abs(tanth);
183 rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
184 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
186 // space-point resolutions
188 sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
189 sigmaZ =AliTPCcluster::SigmaZ2(rpad,tanth );
193 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
194 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
196 // temporary protection
198 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
199 if(sigmaZ < 0.) sigmaZ=0.4e-3;
200 if(clRphi < 0.) clRphi=2.5e-3;
201 if(clZ < 0.) clZ=2.5e-5;
206 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
207 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
209 Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
210 Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
211 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
212 Float_t alpha=(isec < param->GetNInnerSector()) ?
213 param->GetInnerAngle() : param->GetOuterAngle();
214 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
215 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
216 xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
217 if (TMath::Abs(xyz[1])>param->GetZLength()) xyz[1]=tpcHit->Z();
218 xyz[2]=sigmaRphi; // fSigmaY2
219 xyz[3]=sigmaZ; // fSigmaZ2
220 xyz[4]=tpcHit->fQ; // q
222 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
223 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
225 Int_t tracks[3]={tpcHit->Track(), -1, -1};
226 AliTPCcluster cluster(tracks,xyz);
228 clrow->InsertCluster(&cluster); nclusters++;
230 tpcHit = (AliTPChit*)tpc->NextHit();
233 } // end of loop over hits
235 } // end of loop over tracks
237 Int_t nrows=param->GetNRow(isec);
238 for (Int_t irow=0; irow<nrows; irow++) {
239 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
240 if (!clrow) continue;
241 carray.StoreRow(isec,irow);
242 carray.ClearRow(isec,irow);
245 } // end of loop over sectors
247 // cerr<<"Number of made clusters : "<<nclusters<<" \n";
248 loader->WriteRecPoints("OVERWRITE");
253 //_________________________________________________________________
254 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
255 AliTPCClustersArray* clustersArray,
258 //--------------------------------------------------------
259 //calculate exact cross point of track and given pad row
260 //resulting values are expressed in "digit" coordinata
261 //--------------------------------------------------------
263 //-----------------------------------------------------------------
264 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
265 //-----------------------------------------------------------------
267 if (clustersArray==0){
270 AliLoader* loader = runLoader->GetLoader("TPCLoader");
272 AliError("No TPC loader found");
275 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
276 AliRun* aliRun = runLoader->GetAliRun();
278 AliError("Couldn't get AliRun object");
281 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
283 AliError("Couldn't get TPC detector");
286 AliTPCParam* param = tpc->GetParam();
288 AliError("No TPC parameters available");
292 TParticle *particle; // pointer to a given particle
293 AliTPChit *tpcHit; // pointer to a sigle TPC hit
294 // Int_t sector,nhits;
296 const Int_t kcmaxhits=30000;
297 TVector * xxxx = new TVector(kcmaxhits*4);
298 TVector & xxx = *xxxx;
299 Int_t maxhits = kcmaxhits;
300 //construct array for each padrow
301 for (Int_t i=0; i<param->GetNRow(isec);i++)
302 clustersArray->CreateRow(isec,i);
304 //---------------------------------------------------------------
305 // Get the access to the tracks
306 //---------------------------------------------------------------
308 TTree *tH = loader->TreeH();
310 AliFatal("Can not find TreeH in folder");
312 tpc->SetTreeAddress();
314 Stat_t ntracks = tH->GetEntries();
315 Int_t npart = aliRun->GetMCApp()->GetNtrack();
318 if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
319 else branch = tH->GetBranch("TPC");
321 //------------------------------------------------------------
323 //------------------------------------------------------------
325 for(Int_t track=0;track<ntracks;track++){
326 Bool_t isInSector=kTRUE;
328 isInSector = tpc->TrackInVolume(isec,track);
329 if (!isInSector) continue;
331 branch->GetEntry(track); // get next track
333 // Get number of the TPC hits and a pointer
337 Int_t currentIndex=0;
338 Int_t lastrow=-1; //last writen row
342 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
345 Int_t sector=tpcHit->fSector; // sector number
347 tpcHit = (AliTPChit*) tpc->NextHit();
351 ipart=tpcHit->Track();
352 if (ipart<npart) particle=aliRun->GetMCApp()->Particle(ipart);
356 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
357 Int_t index[3]={1,isec,0};
358 Int_t currentrow = param->GetPadRow(x,index) ;
359 if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
360 if (lastrow<0) lastrow=currentrow;
361 if (currentrow==lastrow){
362 if ( currentIndex>=maxhits){
364 xxx.ResizeTo(4*maxhits);
366 xxx(currentIndex*4)=x[0];
367 xxx(currentIndex*4+1)=x[1];
368 xxx(currentIndex*4+2)=x[2];
369 xxx(currentIndex*4+3)=tpcHit->fQ;
385 for (Int_t index=0;index<currentIndex;index++){
387 x=x2=x3=x4=xxx(index*4);
395 sumy+=xxx(index*4+1);
396 sumxy+=xxx(index*4+1)*x;
397 sumx2y+=xxx(index*4+1)*x2;
398 sumz+=xxx(index*4+2);
399 sumxz+=xxx(index*4+2)*x;
400 sumx2z+=xxx(index*4+2)*x2;
401 sumq+=xxx(index*4+3);
403 Float_t centralPad = (param->GetNPads(isec,lastrow)-1)/2;
404 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
405 sumx2*(sumx*sumx3-sumx2*sumx2);
407 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
408 sumx2*(sumxy*sumx3-sumx2y*sumx2);
409 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
410 sumx2*(sumxz*sumx3-sumx2z*sumx2);
412 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
413 sumx2*(sumx*sumx2y-sumx2*sumxy);
414 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
415 sumx2*(sumx*sumx2z-sumx2*sumxz);
417 if (TMath::Abs(det)<0.00001){
418 tpcHit = (AliTPChit*)tpc->NextHit();
422 Float_t y=detay/det+centralPad;
424 Float_t by=detby/det; //y angle
425 Float_t bz=detbz/det; //z angle
426 sumy/=Float_t(currentIndex);
427 sumz/=Float_t(currentIndex);
429 AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
431 AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
432 // AliComplexCluster cl;
438 cl->fTracks[0]=ipart;
442 } //end of calculating cluster for given row
445 tpcHit = (AliTPChit*)tpc->NextHit();
446 } // end of loop over hits
447 } // end of loop over tracks
448 //write padrows to tree
449 for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
450 clustersArray->StoreRow(isec,ii);
451 clustersArray->ClearRow(isec,ii);
458 void AliTPCFast::FindTrackHitsIntersection(AliRunLoader* runLoader,
459 AliTPCClustersArray* clustersArray) const
463 //fill clones array with intersection of current point with the
465 AliLoader* loader = runLoader->GetLoader("TPCLoader");
467 AliError("No TPC loader found");
470 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
471 AliRun* aliRun = runLoader->GetAliRun();
473 AliError("Couldn't get AliRun object");
476 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
478 AliError("Couldn't get TPC detector");
481 AliTPCParam* param = tpc->GetParam();
483 AliError("No TPC parameters available");
490 const Int_t kcmaxhits=30000;
491 TVector * xxxx = new TVector(kcmaxhits*4);
492 TVector & xxx = *xxxx;
493 Int_t maxhits = kcmaxhits;
496 AliTPChit * tpcHit=0;
497 tpcHit = (AliTPChit*)tpc->FirstHit2(-1);
498 Int_t currentIndex=0;
499 Int_t lastrow=-1; //last writen row
502 if (tpcHit==0) continue;
503 sector=tpcHit->fSector; // sector number
504 ipart=tpcHit->Track();
508 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
509 Int_t index[3]={1,sector,0};
510 Int_t currentrow = param->GetPadRow(x,index) ;
511 if (currentrow<0) continue;
512 if (lastrow<0) lastrow=currentrow;
513 if (currentrow==lastrow){
514 if ( currentIndex>=maxhits){
516 xxx.ResizeTo(4*maxhits);
518 xxx(currentIndex*4)=x[0];
519 xxx(currentIndex*4+1)=x[1];
520 xxx(currentIndex*4+2)=x[2];
521 xxx(currentIndex*4+3)=tpcHit->fQ;
537 for (Int_t index=0;index<currentIndex;index++){
539 x=x2=x3=x4=xxx(index*4);
547 sumy+=xxx(index*4+1);
548 sumxy+=xxx(index*4+1)*x;
549 sumx2y+=xxx(index*4+1)*x2;
550 sumz+=xxx(index*4+2);
551 sumxz+=xxx(index*4+2)*x;
552 sumx2z+=xxx(index*4+2)*x2;
553 sumq+=xxx(index*4+3);
555 Float_t centralPad = (param->GetNPads(sector,lastrow)-1)/2;
556 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
557 sumx2*(sumx*sumx3-sumx2*sumx2);
559 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
560 sumx2*(sumxy*sumx3-sumx2y*sumx2);
561 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
562 sumx2*(sumxz*sumx3-sumx2z*sumx2);
564 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
565 sumx2*(sumx*sumx2y-sumx2*sumxy);
566 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
567 sumx2*(sumx*sumx2z-sumx2*sumxz);
569 Float_t y=detay/det+centralPad;
571 Float_t by=detby/det; //y angle
572 Float_t bz=detbz/det; //z angle
573 sumy/=Float_t(currentIndex);
574 sumz/=Float_t(currentIndex);
576 AliComplexCluster cl;
584 AliTPCClustersRow * row = (clustersArray->GetRow(sector,lastrow));
585 if (row!=0) row->InsertCluster(&cl);
588 } //end of calculating cluster for given row
590 } // end of loop over hits