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>
29 #include "AliRunLoader.h"
30 #include "AliLoader.h"
34 #include "AliTPCParam.h"
35 #include "AliTPCClustersArray.h"
36 #include "AliTPCClustersRow.h"
37 #include "AliTPCcluster.h"
38 #include "AliComplexCluster.h"
39 #include "AliTPCFast.h"
43 //____________________________________________________________________
44 AliTPCFast::AliTPCFast(const AliTPCFast ¶m)
45 :TObject(param),fParam(0)
49 // copy constructor - dummy
51 fParam = param.fParam;
53 AliTPCFast & AliTPCFast::operator =(const AliTPCFast & param)
56 // assignment operator - dummy
58 if (this == ¶m) return (*this);
64 //_____________________________________________________________________________
65 void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
67 //--------------------------------------------------------
68 // TPC simple cluster generator from hits
69 // obtained from the TPC Fast Simulator
70 // The point errors are taken from the parametrization
71 //--------------------------------------------------------
73 //-----------------------------------------------------------------
74 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
75 //-----------------------------------------------------------------
76 // Adopted to Marian's cluster data structure by I.Belikov, CERN,
77 // Jouri.Belikov@cern.ch
78 //----------------------------------------------------------------
80 /////////////////////////////////////////////////////////////////////////////
82 //---------------------------------------------------------------------
83 // ALICE TPC Cluster Parameters
84 //--------------------------------------------------------------------
88 // Cluster width in rphi
89 const Float_t kACrphi=0.18322;
90 const Float_t kBCrphi=0.59551e-3;
91 const Float_t kCCrphi=0.60952e-1;
93 const Float_t kACz=0.19081;
94 const Float_t kBCz=0.55938e-3;
95 const Float_t kCCz=0.30428;
98 AliLoader* loader = runLoader->GetLoader("TPCLoader");
100 AliError("No TPC loader found");
103 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
104 AliRun* aliRun = runLoader->GetAliRun();
106 AliError("Couldn't get AliRun object");
109 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
111 AliError("Couldn't get TPC detector");
114 AliTPCParam* param = tpc->GetParam();
116 AliError("No TPC parameters available");
120 //if(fDefaults == 0) SetDefaults();
122 Float_t sigmaRphi,sigmaZ,clRphi,clZ;
124 TParticle *particle; // pointer to a given particle
125 AliTPChit *tpcHit; // pointer to a sigle TPC hit
129 Float_t pl,pt,tanth,rpad,ratio;
132 //---------------------------------------------------------------
133 // Get the access to the tracks
134 //---------------------------------------------------------------
136 TTree *tH = loader->TreeH();
138 AliFatal("Can not find TreeH in folder");
141 tpc->SetTreeAddress();
143 Stat_t ntracks = tH->GetEntries();
145 //Switch to the output file
147 if (loader->TreeR() == 0x0) loader->MakeTree("R");
149 AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
151 runLoader->CdGAFile();
152 //param->Write(param->GetTitle());
154 AliTPCClustersArray carray;
156 carray.SetClusterType("AliTPCcluster");
157 carray.MakeTree(loader->TreeR());
159 Int_t nclusters=0; //cluster counter
161 //------------------------------------------------------------
162 // Loop over all sectors (72 sectors for 20 deg
163 // segmentation for both lower and upper sectors)
164 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
165 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
167 // First cluster for sector 0 starts at "0"
168 //------------------------------------------------------------
170 for(Int_t isec=0;isec<param->GetNSector();isec++){
172 param->AdjustCosSin(isec,cph,sph);
174 //------------------------------------------------------------
176 //------------------------------------------------------------
178 for(Int_t track=0;track<ntracks;track++){
182 // Get number of the TPC hits
184 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
190 if (tpcHit->fQ == 0.) {
191 tpcHit = (AliTPChit*) tpc->NextHit();
192 continue; //information about track (I.Belikov)
194 sector=tpcHit->fSector; // sector number
197 tpcHit = (AliTPChit*) tpc->NextHit();
200 ipart=tpcHit->Track();
201 particle=aliRun->GetMCApp()->Particle(ipart);
204 if(pt < 1.e-9) pt=1.e-9;
206 tanth = TMath::Abs(tanth);
207 rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
208 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
210 // space-point resolutions
212 sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
213 sigmaZ =AliTPCcluster::SigmaZ2(rpad,tanth );
217 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
218 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
220 // temporary protection
222 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
223 if(sigmaZ < 0.) sigmaZ=0.4e-3;
224 if(clRphi < 0.) clRphi=2.5e-3;
225 if(clZ < 0.) clZ=2.5e-5;
230 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
231 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
233 Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
234 Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
235 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
236 Float_t alpha=(isec < param->GetNInnerSector()) ?
237 param->GetInnerAngle() : param->GetOuterAngle();
238 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
239 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
240 xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
241 if (TMath::Abs(xyz[1])>param->GetZLength(isec)) xyz[1]=tpcHit->Z();
242 xyz[2]=sigmaRphi; // fSigmaY2
243 xyz[3]=sigmaZ; // fSigmaZ2
244 xyz[4]=tpcHit->fQ; // q
246 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
247 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
249 Int_t tracks[3]={tpcHit->Track(), -1, -1};
250 AliTPCcluster cluster(tracks,xyz);
252 clrow->InsertCluster(&cluster); nclusters++;
254 tpcHit = (AliTPChit*)tpc->NextHit();
257 } // end of loop over hits
259 } // end of loop over tracks
261 Int_t nrows=param->GetNRow(isec);
262 for (Int_t irow=0; irow<nrows; irow++) {
263 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
264 if (!clrow) continue;
265 carray.StoreRow(isec,irow);
266 carray.ClearRow(isec,irow);
269 } // end of loop over sectors
271 // cerr<<"Number of made clusters : "<<nclusters<<" \n";
272 loader->WriteRecPoints("OVERWRITE");
279 //_________________________________________________________________
283 void AliTPCFast::Hits2ExactClusters(AliRunLoader* runLoader) const{
287 AliLoader* loader = runLoader->GetLoader("TPCLoader");
289 AliError("No TPC loader found");
292 AliTPC* tpc = (AliTPC*) gAlice->GetDetector("TPC");
294 AliError("Couldn't get TPC detector");
297 AliTPCParam* param = tpc->GetParam();
299 AliError("No TPC parameters available");
302 //---------------------------------------------------------------
303 // Get the access to the tracks
304 //---------------------------------------------------------------
306 TTree *tH = loader->TreeH();
307 if (tH == 0x0) { runLoader->LoadHits("TPC"); loader->TreeH();}
309 tpc->SetTreeAddress();
312 //Switch to the output file
314 if (loader->TreeR() == 0x0) loader->MakeTree("R");
316 AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
318 runLoader->CdGAFile();
319 //param->Write(param->GetTitle());
321 AliTPCClustersArray carray;
323 carray.SetClusterType("AliTPCclusterMI");
324 carray.MakeTree(loader->TreeR());
326 //------------------------------------------------------------
327 // Loop over all sectors (72 sectors for 20 deg)
328 //------------------------------------------------------------
330 for(Int_t isec=0;isec<param->GetNSector();isec++){
331 Hits2ExactClustersSector(runLoader, &carray, isec);
332 } // end of loop over sectors
333 loader->WriteRecPoints("OVERWRITE");
339 //_________________________________________________________________
340 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
341 AliTPCClustersArray* clustersArray,
344 //--------------------------------------------------------
345 //calculate exact cross point of track and given pad row
346 //resulting values are expressed in "local" coordinata
347 //the sigmaY2 and sigmaZ2 of cluster are the shape of cluster parameters
348 // - thay are used later on for error parameterization
349 //--------------------------------------------------------
351 //-----------------------------------------------------------------
352 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
353 //-----------------------------------------------------------------
355 if (clustersArray==0){
358 AliLoader* loader = runLoader->GetLoader("TPCLoader");
360 AliError("No TPC loader found");
363 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
364 AliRun* aliRun = runLoader->GetAliRun();
366 AliError("Couldn't get AliRun object");
369 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
371 AliError("Couldn't get TPC detector");
374 AliTPCParam* param = tpc->GetParam();
376 AliError("No TPC parameters available");
381 AliTPChit *tpcHit; // pointer to a sigle TPC hit
382 // Int_t sector,nhits;
384 const Int_t kcmaxhits=30000;
385 TVector * xxxx = new TVector(kcmaxhits*4);
386 TVector & xxx = *xxxx;
387 Int_t maxhits = kcmaxhits;
388 //construct array for each padrow
389 for (Int_t i=0; i<param->GetNRow(isec);i++)
390 clustersArray->CreateRow(isec,i);
392 //---------------------------------------------------------------
393 // Get the access to the tracks
394 //---------------------------------------------------------------
396 TTree *tH = loader->TreeH();
398 AliFatal("Can not find TreeH in folder");
401 tpc->SetTreeAddress();
403 Stat_t ntracks = tH->GetEntries();
406 if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
407 else branch = tH->GetBranch("TPC");
409 //------------------------------------------------------------
411 //------------------------------------------------------------
413 for(Int_t track=0;track<ntracks;track++){
414 Bool_t isInSector=kTRUE;
416 isInSector = tpc->TrackInVolume(isec,track);
417 if (!isInSector) continue;
419 branch->GetEntry(track); // get next track
423 Int_t currentIndex=0;
424 Int_t lastrow=-1; //last writen row
428 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
431 Int_t sector=tpcHit->fSector; // sector number
433 tpcHit = (AliTPChit*) tpc->NextHit();
437 ipart=tpcHit->Track();
441 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
442 Int_t index[3]={1,isec,0};
443 Int_t currentrow = param->GetPadRow(x,index);
445 if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
446 if (lastrow<0) lastrow=currentrow;
449 if (currentrow!=lastrow){
462 for (Int_t cindex=0;cindex<currentIndex;cindex++){
464 x1=x2=x3=x4=xxx(cindex*4);
472 sumy+=xxx(cindex*4+1);
473 sumxy+=xxx(cindex*4+1)*x1;
474 sumx2y+=xxx(cindex*4+1)*x2;
475 sumz+=xxx(cindex*4+2);
476 sumxz+=xxx(cindex*4+2)*x1;
477 sumx2z+=xxx(cindex*4+2)*x2;
478 sumq+=xxx(cindex*4+3);
480 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
481 sumx2*(sumx*sumx3-sumx2*sumx2);
483 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
484 sumx2*(sumxy*sumx3-sumx2y*sumx2);
485 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
486 sumx2*(sumxz*sumx3-sumx2z*sumx2);
488 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
489 sumx2*(sumx*sumx2y-sumx2*sumxy);
490 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
491 sumx2*(sumx*sumx2z-sumx2*sumxz);
493 if (TMath::Abs(det)<0.00000000000000001){
494 tpcHit = (AliTPChit*)tpc->NextHit();
500 Float_t by=detby/det; //y angle
501 Float_t bz=detbz/det; //z angle
502 sumy/=Float_t(currentIndex);
503 sumz/=Float_t(currentIndex);
508 Float_t sign = (tpcHit->Z()>0)? 1.:-1.;
509 y = (y+0.5)*param->GetPadPitchWidth(isec);
510 z = z*param->GetZWidth();
513 Double_t sigmay2 = z*fParam->GetDiffL();
515 Float_t angulary = by*param->GetPadPitchLength(isec,lastrow);
518 sigmay2 +=angulary+0.25*param->GetPadPitchWidth(isec)*param->GetPadPitchWidth(isec);
521 Double_t sigmaz2 = z*fParam->GetDiffT();
523 Float_t angularz = bz*param->GetPadPitchLength(isec,lastrow);
526 sigmaz2 +=angularz+0.25;
528 sigmaz2 = TMath::Min(sigmaz2,1.);
529 sigmay2 = TMath::Min(sigmay2,1.);
532 z = sign*(param->GetZLength(isec) - z);
533 if (TMath::Abs(z)< param->GetZLength(isec)-1){
534 AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
536 AliTPCclusterMI* cl = new((AliTPCclusterMI*)row->Append()) AliTPCclusterMI ;
537 // AliTPCclusterMI cl;
541 if (TMath::Abs(sumq)<1000){
542 cl->SetMax(Short_t(sumq));
547 cl->SetSigmaZ2(sigmaz2);
548 cl->SetSigmaY2(sigmay2);
549 cl->SetLabel(ipart,0);
555 } //end of calculating cluster for given row
558 } // end of crossing rows
560 if ( currentIndex>=maxhits){
562 xxx.ResizeTo(4*maxhits);
564 xxx(currentIndex*4)=x[0];
565 xxx(currentIndex*4+1)=x[1];
566 xxx(currentIndex*4+2)=x[2];
567 xxx(currentIndex*4+3)=tpcHit->fQ;
569 tpcHit = (AliTPChit*)tpc->NextHit();
570 } // end of loop over hits
571 } // end of loop over tracks
572 //write padrows to tree
573 for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
574 clustersArray->StoreRow(isec,ii);
575 clustersArray->ClearRow(isec,ii);