7e074f739758bcf9bbbff336d6ac033abbbe4630
[u/mrichter/AliRoot.git] / TPC / AliTPCFast.cxx
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>
25 #include <TVector.h>
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"
37 #include "AliLog.h"
38
39 ClassImp(AliTPCFast)
40
41
42 //_____________________________________________________________________________
43 void 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) {
78     AliError("No TPC loader found");
79     return;
80   }
81   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
82   AliRun* aliRun = runLoader->GetAliRun();
83   if (!aliRun) {
84     AliError("Couldn't get AliRun object");
85     return;
86   }
87   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
88   if (!tpc) {
89     AliError("Couldn't get TPC detector");
90     return;
91   }
92   AliTPCParam* param = tpc->GetParam();
93   if (!param) {
94     AliError("No TPC parameters available");
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)
116     AliFatal("Can not find TreeH in folder");
117
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   
126   AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
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){
166  
167          if (tpcHit->fQ == 0.) {
168            tpcHit = (AliTPChit*) tpc->NextHit();
169            continue; //information about track (I.Belikov)
170          }
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         
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
254 //_________________________________________________________________
255 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
256                                           AliTPCClustersArray* clustersArray,
257                                           Int_t isec) const
258 {
259   //--------------------------------------------------------
260   //calculate exact cross point of track and given pad row
261   //resulting values are expressed in "digit" coordinata
262   //--------------------------------------------------------
263
264   //-----------------------------------------------------------------
265   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
266   //-----------------------------------------------------------------
267   //
268   if (clustersArray==0){    
269     return;
270   }
271   AliLoader* loader = runLoader->GetLoader("TPCLoader");
272   if (!loader) {
273     AliError("No TPC loader found");
274     return;
275   }
276   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
277   AliRun* aliRun = runLoader->GetAliRun();
278   if (!aliRun) {
279     AliError("Couldn't get AliRun object");
280     return;
281   }
282   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
283   if (!tpc) {
284     AliError("Couldn't get TPC detector");
285     return;
286   }
287   AliTPCParam* param = tpc->GetParam();
288   if (!param) {
289     AliError("No TPC parameters available");
290     return;
291   }
292   //
293   TParticle *particle; // pointer to a given particle
294   AliTPChit *tpcHit; // pointer to a sigle TPC hit
295   //  Int_t sector,nhits;
296   Int_t ipart;
297   const Int_t kcmaxhits=30000;
298   TVector * xxxx = new TVector(kcmaxhits*4);
299   TVector & xxx = *xxxx;
300   Int_t maxhits = kcmaxhits;
301   //construct array for each padrow
302   for (Int_t i=0; i<param->GetNRow(isec);i++) 
303     clustersArray->CreateRow(isec,i);
304   
305   //---------------------------------------------------------------
306   //  Get the access to the tracks 
307   //---------------------------------------------------------------
308   
309   TTree *tH = loader->TreeH();
310   if (tH == 0x0)
311     AliFatal("Can not find TreeH in folder");
312
313   tpc->SetTreeAddress();
314
315   Stat_t ntracks = tH->GetEntries();
316   Int_t npart = aliRun->GetMCApp()->GetNtrack();
317   //MI change
318   TBranch * branch=0;
319   if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
320   else branch = tH->GetBranch("TPC");
321
322   //------------------------------------------------------------
323   // Loop over tracks
324   //------------------------------------------------------------
325
326   for(Int_t track=0;track<ntracks;track++){ 
327     Bool_t isInSector=kTRUE;
328     tpc->ResetHits();
329      isInSector = tpc->TrackInVolume(isec,track);
330     if (!isInSector) continue;
331     //MI change
332     branch->GetEntry(track); // get next track
333     //
334     //  Get number of the TPC hits and a pointer
335     //  to the particles
336     // Loop over hits
337     //
338     Int_t currentIndex=0;
339     Int_t lastrow=-1;  //last writen row
340
341     //M.I. changes
342
343     tpcHit = (AliTPChit*)tpc->FirstHit(-1);
344     while(tpcHit){
345       
346       Int_t sector=tpcHit->fSector; // sector number
347       if(sector != isec){
348         tpcHit = (AliTPChit*) tpc->NextHit();
349         continue; 
350       }
351
352       ipart=tpcHit->Track();
353       if (ipart<npart) particle=aliRun->GetMCApp()->Particle(ipart);
354       
355       //find row number
356
357       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
358       Int_t    index[3]={1,isec,0};
359       Int_t    currentrow = param->GetPadRow(x,index) ; 
360       if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
361       if (lastrow<0) lastrow=currentrow;
362       if (currentrow==lastrow){
363         if ( currentIndex>=maxhits){
364           maxhits+=kcmaxhits;
365           xxx.ResizeTo(4*maxhits);
366         }     
367         xxx(currentIndex*4)=x[0];
368         xxx(currentIndex*4+1)=x[1];
369         xxx(currentIndex*4+2)=x[2];     
370         xxx(currentIndex*4+3)=tpcHit->fQ;
371         currentIndex++; 
372       }
373       else 
374         if (currentIndex>2){
375           Float_t sumx=0;
376           Float_t sumx2=0;
377           Float_t sumx3=0;
378           Float_t sumx4=0;
379           Float_t sumy=0;
380           Float_t sumxy=0;
381           Float_t sumx2y=0;
382           Float_t sumz=0;
383           Float_t sumxz=0;
384           Float_t sumx2z=0;
385           Float_t sumq=0;
386           for (Int_t index=0;index<currentIndex;index++){
387             Float_t x,x2,x3,x4;
388             x=x2=x3=x4=xxx(index*4);
389             x2*=x;
390             x3*=x2;
391             x4*=x3;
392             sumx+=x;
393             sumx2+=x2;
394             sumx3+=x3;
395             sumx4+=x4;
396             sumy+=xxx(index*4+1);
397             sumxy+=xxx(index*4+1)*x;
398             sumx2y+=xxx(index*4+1)*x2;
399             sumz+=xxx(index*4+2);
400             sumxz+=xxx(index*4+2)*x;
401             sumx2z+=xxx(index*4+2)*x2;   
402             sumq+=xxx(index*4+3);
403           }
404           Float_t centralPad = (param->GetNPads(isec,lastrow)-1)/2;
405           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
406             sumx2*(sumx*sumx3-sumx2*sumx2);
407           
408           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
409             sumx2*(sumxy*sumx3-sumx2y*sumx2);
410           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
411             sumx2*(sumxz*sumx3-sumx2z*sumx2);
412           
413           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
414             sumx2*(sumx*sumx2y-sumx2*sumxy);
415           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
416             sumx2*(sumx*sumx2z-sumx2*sumxz);
417           
418           if (TMath::Abs(det)<0.00001){
419              tpcHit = (AliTPChit*)tpc->NextHit();
420             continue;
421           }
422         
423           Float_t y=detay/det+centralPad;
424           Float_t z=detaz/det;  
425           Float_t by=detby/det; //y angle
426           Float_t bz=detbz/det; //z angle
427           sumy/=Float_t(currentIndex);
428           sumz/=Float_t(currentIndex);
429
430           AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
431           if (row!=0) {
432             AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
433             //    AliComplexCluster cl;
434             cl->fX=z;
435             cl->fY=y;
436             cl->fQ=sumq;
437             cl->fSigmaX2=bz;
438             cl->fSigmaY2=by;
439             cl->fTracks[0]=ipart;
440           }
441           currentIndex=0;
442           lastrow=currentrow;
443         } //end of calculating cluster for given row
444         
445         
446       tpcHit = (AliTPChit*)tpc->NextHit();
447     } // end of loop over hits
448   }   // end of loop over tracks 
449   //write padrows to tree 
450   for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
451     clustersArray->StoreRow(isec,ii);    
452     clustersArray->ClearRow(isec,ii);        
453   }
454   xxxx->Delete();
455  
456 }
457
458
459 void AliTPCFast::FindTrackHitsIntersection(AliRunLoader* runLoader,
460                                            AliTPCClustersArray* clustersArray) const
461 {
462
463   //
464   //fill clones array with intersection of current point with the
465   //middle of the row
466   AliLoader* loader = runLoader->GetLoader("TPCLoader");
467   if (!loader) {
468     AliError("No TPC loader found");
469     return;
470   }
471   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
472   AliRun* aliRun = runLoader->GetAliRun();
473   if (!aliRun) {
474     AliError("Couldn't get AliRun object");
475     return;
476   }
477   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
478   if (!tpc) {
479     AliError("Couldn't get TPC detector");
480     return;
481   }
482   AliTPCParam* param = tpc->GetParam();
483   if (!param) {
484     AliError("No TPC parameters available");
485     return;
486   }
487
488   Int_t sector;
489   Int_t ipart;
490   
491   const Int_t kcmaxhits=30000;
492   TVector * xxxx = new TVector(kcmaxhits*4);
493   TVector & xxx = *xxxx;
494   Int_t maxhits = kcmaxhits;
495       
496   //
497   AliTPChit * tpcHit=0;
498   tpcHit = (AliTPChit*)tpc->FirstHit2(-1);
499   Int_t currentIndex=0;
500   Int_t lastrow=-1;  //last writen row
501
502   while (tpcHit){
503     if (tpcHit==0) continue;
504     sector=tpcHit->fSector; // sector number
505     ipart=tpcHit->Track();
506     
507     //find row number
508     
509     Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
510     Int_t    index[3]={1,sector,0};
511     Int_t    currentrow = param->GetPadRow(x,index) ;   
512     if (currentrow<0) continue;
513     if (lastrow<0) lastrow=currentrow;
514     if (currentrow==lastrow){
515       if ( currentIndex>=maxhits){
516         maxhits+=kcmaxhits;
517         xxx.ResizeTo(4*maxhits);
518       }     
519       xxx(currentIndex*4)=x[0];
520       xxx(currentIndex*4+1)=x[1];
521       xxx(currentIndex*4+2)=x[2];       
522       xxx(currentIndex*4+3)=tpcHit->fQ;
523       currentIndex++;   
524     }
525     else 
526       if (currentIndex>2){
527         Float_t sumx=0;
528         Float_t sumx2=0;
529         Float_t sumx3=0;
530         Float_t sumx4=0;
531         Float_t sumy=0;
532         Float_t sumxy=0;
533         Float_t sumx2y=0;
534         Float_t sumz=0;
535         Float_t sumxz=0;
536         Float_t sumx2z=0;
537         Float_t sumq=0;
538         for (Int_t index=0;index<currentIndex;index++){
539           Float_t x,x2,x3,x4;
540           x=x2=x3=x4=xxx(index*4);
541           x2*=x;
542           x3*=x2;
543           x4*=x3;
544           sumx+=x;
545           sumx2+=x2;
546           sumx3+=x3;
547           sumx4+=x4;
548           sumy+=xxx(index*4+1);
549           sumxy+=xxx(index*4+1)*x;
550           sumx2y+=xxx(index*4+1)*x2;
551           sumz+=xxx(index*4+2);
552           sumxz+=xxx(index*4+2)*x;
553           sumx2z+=xxx(index*4+2)*x2;     
554           sumq+=xxx(index*4+3);
555         }
556         Float_t centralPad = (param->GetNPads(sector,lastrow)-1)/2;
557         Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
558           sumx2*(sumx*sumx3-sumx2*sumx2);
559         
560         Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
561           sumx2*(sumxy*sumx3-sumx2y*sumx2);
562         Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
563           sumx2*(sumxz*sumx3-sumx2z*sumx2);
564         
565         Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
566           sumx2*(sumx*sumx2y-sumx2*sumxy);
567         Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
568           sumx2*(sumx*sumx2z-sumx2*sumxz);
569         
570         Float_t y=detay/det+centralPad;
571         Float_t z=detaz/det;    
572         Float_t by=detby/det; //y angle
573         Float_t bz=detbz/det; //z angle
574         sumy/=Float_t(currentIndex);
575         sumz/=Float_t(currentIndex);
576         
577         AliComplexCluster cl;
578         cl.fX=z;
579         cl.fY=y;
580         cl.fQ=sumq;
581         cl.fSigmaX2=bz;
582         cl.fSigmaY2=by;
583         cl.fTracks[0]=ipart;
584         
585         AliTPCClustersRow * row = (clustersArray->GetRow(sector,lastrow));
586         if (row!=0) row->InsertCluster(&cl);
587         currentIndex=0;
588         lastrow=currentrow;
589       } //end of calculating cluster for given row
590                 
591   } // end of loop over hits
592   xxxx->Delete();
593
594 }
595