ATO-98 - function returns also status value
[u/mrichter/AliRoot.git] / PMD / AliPMDtracker.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 //-----------------------------------------------------//
17 //                                                     //
18 //           Date   : March 25 2004                    //
19 //  This reads the file PMD.RecPoints.root(TreeR),     //
20 //  calls the Clustering algorithm and stores the      //
21 //  clustering output in PMD.RecPoints.root(TreeR)     // 
22 //                                                     //
23 //-----------------------------------------------------//
24
25 #include <Riostream.h>
26 #include <TMath.h>
27 #include <TTree.h>
28 #include <TObjArray.h>
29 #include <TClonesArray.h>
30 #include <TFile.h>
31 #include <TBranch.h>
32 #include <TNtuple.h>
33 #include <TParticle.h>
34
35 #include <TGeoMatrix.h>
36
37 #include "AliGeomManager.h"
38
39 #include "AliPMDcluster.h"
40 #include "AliPMDclupid.h"
41 #include "AliPMDrecpoint1.h"
42 #include "AliPMDrecdata.h"
43 #include "AliPMDrechit.h"
44 #include "AliPMDUtility.h"
45 #include "AliPMDDiscriminator.h"
46 #include "AliPMDEmpDiscriminator.h"
47 #include "AliPMDtracker.h"
48
49 #include "AliESDPmdTrack.h"
50 #include "AliESDEvent.h"
51 #include "AliLog.h"
52
53 ClassImp(AliPMDtracker)
54
55 AliPMDtracker::AliPMDtracker():
56   fTreeR(0),
57   fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
58   fRechits(new TClonesArray("AliPMDrechit", 10)),
59   fPMDcontin(new TObjArray()),
60   fPMDcontout(new TObjArray()),
61   fPMDutil(new AliPMDUtility()),
62   fPMDrecpoint(0),
63   fPMDclin(0),
64   fPMDclout(0),
65   fXvertex(0.),
66   fYvertex(0.),
67   fZvertex(0.),
68   fSigmaX(0.),
69   fSigmaY(0.),
70   fSigmaZ(0.)
71 {
72   //
73   // Default Constructor
74   //
75 }
76 //--------------------------------------------------------------------//
77 AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
78   TObject(/* tracker */),
79   fTreeR(0),
80   fRecpoints(NULL),
81   fRechits(NULL),
82   fPMDcontin(NULL),
83   fPMDcontout(NULL),
84   fPMDutil(NULL),
85   fPMDrecpoint(0),
86   fPMDclin(0),
87   fPMDclout(0),
88   fXvertex(0.),
89   fYvertex(0.),
90   fZvertex(0.),
91   fSigmaX(0.),
92   fSigmaY(0.),
93   fSigmaZ(0.)
94 {
95   // copy constructor
96   AliError("Copy constructor not allowed");
97 }
98
99 //--------------------------------------------------------------------//
100 AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & /* tracker */)
101 {
102  // assignment operator
103   AliError("Assignment operator not allowed");
104   return *this;
105 }
106
107 //--------------------------------------------------------------------//
108 AliPMDtracker::~AliPMDtracker()
109 {
110   // Destructor
111   if (fRecpoints)
112     {
113       fRecpoints->Clear();
114     }
115   if (fRechits)
116     {
117       fRechits->Clear();
118     }
119
120   if (fPMDcontin)
121     {
122       fPMDcontin->Delete();
123       delete fPMDcontin;
124       fPMDcontin=0;
125       
126     }
127   if (fPMDcontout)
128   {
129       fPMDcontout->Delete();
130       delete fPMDcontout;
131       fPMDcontout=0;
132
133     }
134   delete fPMDutil;
135 }
136 //--------------------------------------------------------------------//
137 void AliPMDtracker::LoadClusters(TTree *treein)
138 {
139   // Load the Reconstructed tree
140   fTreeR = treein;
141 }
142 //--------------------------------------------------------------------//
143 void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
144 {
145   // Converts digits to recpoints after running clustering
146   // algorithm on CPV plane and PREshower plane
147   //
148
149   Int_t   idet = 0;
150   Int_t   ismn = 0;
151   Int_t   trackno = 1, trackpid = 0;
152   Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
153   
154   Int_t *irow;
155   Int_t *icol;
156   Int_t *itra;
157   Int_t *ipid;
158   Float_t *cadc;
159
160   AliPMDrechit *rechit = 0x0;
161
162   TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
163   if (!branch)
164     {
165       AliError("PMDRecpoint branch not found");
166       return;
167     }
168   branch->SetAddress(&fRecpoints);  
169
170   TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
171   if (!branch1)
172     {
173       AliError("PMDRechit branch not found");
174       return;
175     }
176   branch1->SetAddress(&fRechits);  
177
178   Int_t ncrhit = 0;
179   Int_t   nmodules = (Int_t) branch->GetEntries();
180   
181   AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
182   for (Int_t imodule = 0; imodule < nmodules; imodule++)
183     {
184       branch->GetEntry(imodule); 
185       Int_t nentries = fRecpoints->GetLast();
186       AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
187                       ,nentries));
188
189       for(Int_t ient = 0; ient < nentries+1; ient++)
190         {
191           fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
192           idet        = fPMDrecpoint->GetDetector();
193           ismn        = fPMDrecpoint->GetSMNumber();
194           clusdata[0] = fPMDrecpoint->GetClusX();
195           clusdata[1] = fPMDrecpoint->GetClusY();
196           clusdata[2] = fPMDrecpoint->GetClusADC();
197           clusdata[3] = fPMDrecpoint->GetClusCells();
198           clusdata[4] = fPMDrecpoint->GetClusSigmaX();
199           clusdata[5] = fPMDrecpoint->GetClusSigmaY();
200
201           if (clusdata[4] >= 0. && clusdata[5] >= 0.)
202             { 
203               // extract the associated cell information
204               branch1->GetEntry(ncrhit); 
205               Int_t nenbr1 = fRechits->GetLast() + 1;
206
207               irow = new Int_t[nenbr1];
208               icol = new Int_t[nenbr1];
209               itra = new Int_t[nenbr1];
210               ipid = new Int_t[nenbr1];
211               cadc = new Float_t[nenbr1];
212               
213               for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
214                 {
215                   irow[ient1] = -99;
216                   icol[ient1] = -99;
217                   itra[ient1] = -99;
218                   ipid[ient1] = -99;
219                   cadc[ient1] = 0.;
220                 }
221               for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
222                 {
223                   rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
224                   //irow[ient1] = rechit->GetCellX();
225                   //icol[ient1] = rechit->GetCellY();
226                   itra[ient1] = rechit->GetCellTrack();
227                   ipid[ient1] = rechit->GetCellPid();
228                   cadc[ient1] = rechit->GetCellAdc();
229                 }
230               if (idet == 0)
231                 {
232                   AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
233                                        trackno, trackpid);
234                 }
235               else if (idet == 1)
236                 {
237                   trackno  = itra[0];
238                   trackpid = ipid[0];
239                 }
240
241               delete [] irow;
242               delete [] icol;
243               delete [] itra;
244               delete [] ipid;
245               delete [] cadc;
246
247               fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
248               fPMDcontin->Add(fPMDclin);
249
250               ncrhit++;
251             }
252         }
253     }
254
255   AliPMDEmpDiscriminator pmddiscriminator;
256   pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
257
258   // alignment implemention
259
260   Double_t sectr[4][3] = { {0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
261   TString snsector="PMD/Sector";
262   TString symname;
263   TGeoHMatrix gpmdor;
264   
265   for(Int_t isector=1; isector<=4; isector++)
266     {
267       symname = snsector;
268       symname += isector;
269       TGeoHMatrix *gpmdal = AliGeomManager::GetMatrix(symname);
270       Double_t *tral = gpmdal->GetTranslation();
271
272       AliGeomManager::GetOrigGlobalMatrix(symname, gpmdor);
273       Double_t *tror = gpmdor.GetTranslation();
274       
275       for(Int_t ixyz=0; ixyz<3; ixyz++)
276         {
277           sectr[isector-1][ixyz] = tral[ixyz] - tror[ixyz];
278         }
279     }
280
281   const Float_t kzpos = 361.5;    // middle of the PMD
282
283   Int_t   ix = -1, iy = -1;
284   Int_t   det = 0, smn = 0, trno = 1, trpid = 0, mstat = 0;
285   Float_t xpos = 0., ypos = 0.;
286   Float_t adc = 0., ncell = 0., radx = 0., rady = 0.;
287   Float_t xglobal = 0., yglobal = 0., zglobal = 0;
288   Float_t pid = 0.;
289
290   fPMDutil->ApplyAlignment(sectr);
291
292   Int_t nentries2 = fPMDcontout->GetEntries();
293   AliDebug(1,Form("Number of clusters coming after discrimination = %d"
294                   ,nentries2));
295   for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
296     {
297       fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
298       
299       det   = fPMDclout->GetDetector();
300       smn   = fPMDclout->GetSMN();
301       trno  = fPMDclout->GetClusTrackNo();
302       trpid = fPMDclout->GetClusTrackPid();
303       mstat = fPMDclout->GetClusMatching();
304       xpos  = fPMDclout->GetClusX();
305       ypos  = fPMDclout->GetClusY();
306       adc   = fPMDclout->GetClusADC();
307       ncell = fPMDclout->GetClusCells();
308       radx  = fPMDclout->GetClusSigmaX();
309       // Here in the variable "rady" we are keeping the row and col
310       // of the single isolated cluster having ncell = 1 for offline
311       // calibration
312       
313       if ((radx > 999. && radx < 1000.) && ncell == 1)
314         {
315           if (smn < 12)
316             {
317               ix = (Int_t) (ypos +0.5);
318               iy = (Int_t) xpos;
319             }
320           else if (smn > 12 && smn < 24)
321             {
322               ix = (Int_t) xpos;
323               iy = (Int_t) (ypos +0.5);
324             }
325           rady = (Float_t) (ix*100 + iy);
326         }
327       else
328         {
329           rady  = fPMDclout->GetClusSigmaY();
330         }
331       pid   = fPMDclout->GetClusPID();
332       
333       //
334       /**********************************************************************
335        *    det   : Detector, 0: PRE & 1:CPV                                *
336        *    smn   : Serial Module Number 0 to 23 for each plane             *
337        *    xpos  : x-position of the cluster                               *
338        *    ypos  : y-position of the cluster                               *
339        *            THESE xpos & ypos are not the true xpos and ypos        *
340        *            for some of the unit modules. They are rotated.         *
341        *    adc   : ADC contained in the cluster                            *
342        *    ncell : Number of cells contained in the cluster                *
343        *    rad   : radius of the cluster (1d fit)                          *
344        **********************************************************************/
345       //
346
347
348       if (det == 0)
349         {
350           zglobal = kzpos + 1.65; // PREshower plane
351         }
352       else if (det == 1)
353         {
354           zglobal = kzpos - 1.65; // CPV plane
355         }
356
357       fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
358
359       // Fill ESD
360
361       AliESDPmdTrack *esdpmdtr = new  AliESDPmdTrack();
362
363       esdpmdtr->SetDetector(det);
364       esdpmdtr->SetSmn(smn);
365       esdpmdtr->SetClusterTrackNo(trno);
366       esdpmdtr->SetClusterTrackPid(trpid);
367       esdpmdtr->SetClusterMatching(mstat);
368       
369       esdpmdtr->SetClusterX(xglobal);
370       esdpmdtr->SetClusterY(yglobal);
371       esdpmdtr->SetClusterZ(zglobal);
372       esdpmdtr->SetClusterADC(adc);
373       esdpmdtr->SetClusterCells(ncell);
374       esdpmdtr->SetClusterPID(pid);
375       esdpmdtr->SetClusterSigmaX(radx);
376       esdpmdtr->SetClusterSigmaY(rady);
377
378       event->AddPmdTrack(esdpmdtr);
379       delete esdpmdtr;
380     }
381
382   fPMDcontin->Delete();
383   fPMDcontout->Delete();
384
385 }
386 //--------------------------------------------------------------------//
387 void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
388                                          Int_t *ipid, Float_t *cadc,
389                                          Int_t &trackno, Int_t &trackpid)
390 {
391   // assign the track number and the corresponding pid to a cluster
392   // split cluster part will be done at the time of calculating eff/pur
393
394   Int_t *phentry = new Int_t [nentry];
395   Int_t *hadentry = new Int_t [nentry];
396   Int_t *trpid    = 0x0;
397   Int_t *sortcoord = 0x0;
398   Float_t *trenergy = 0x0;
399
400   Int_t ngtrack = 0;
401   Int_t nhtrack = 0;
402   for (Int_t i = 0; i < nentry; i++)
403     {
404       phentry[i] = -1;
405       hadentry[i] = -1;
406
407       if (ipid[i] == 22)
408         {
409           phentry[ngtrack] = i;
410           ngtrack++;
411         }
412       else if (ipid[i] != 22)
413         {
414           hadentry[nhtrack] = i;
415           nhtrack++;
416         }
417     }
418   
419   Int_t nghadtrack = ngtrack + nhtrack;
420
421   if (ngtrack == 0)
422     {
423       // hadron track
424       // no need of track number, set to -1
425       trackpid = 8;
426       trackno  = -1;
427     }
428   else if (ngtrack >= 1)
429     {
430       // one or more than one photon track + charged track
431       // find out which track deposits maximum energy and
432       // assign that track number and track pid
433
434       trenergy  = new Float_t [nghadtrack];
435       trpid     = new Int_t [nghadtrack];
436       sortcoord = new Int_t [2*nghadtrack];
437       for (Int_t i = 0; i < ngtrack; i++)
438         {
439           trenergy[i] = 0.;
440           trpid[i]    = -1;
441           for (Int_t j = 0; j < nentry; j++)
442             {
443               if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
444                 {
445                   trenergy[i] += cadc[j];
446                   trpid[i]     = 22;
447                 }
448             }
449         }
450       for (Int_t i = ngtrack; i < nghadtrack; i++)
451         {
452           trenergy[i] = 0.;
453           trpid[i]    = -1;
454           for (Int_t j = 0; j < nentry; j++)
455             {
456               if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
457                 {
458                   trenergy[i] += cadc[j];
459                   trpid[i]     = ipid[j];
460                 }
461             }
462         }
463       
464       Bool_t jsort = true;
465       TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
466       
467       Int_t gtr = sortcoord[0];   
468       if (trpid[gtr] == 22)
469         {
470           trackpid = 22;
471           trackno  = itra[phentry[gtr]];   // highest adc track
472         }
473       else
474         {
475           trackpid = 8;
476           trackno = -1;
477         }
478       
479       delete [] trenergy;
480       delete [] trpid;
481       delete [] sortcoord;
482       
483     }   // end of ngtrack >= 1
484
485   delete [] phentry;
486   delete [] hadentry;
487   
488 }
489 //--------------------------------------------------------------------//
490 void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
491 {
492   fXvertex = vtx[0];
493   fYvertex = vtx[1];
494   fZvertex = vtx[2];
495   fSigmaX  = evtx[0];
496   fSigmaY  = evtx[1];
497   fSigmaZ  = evtx[2];
498 }
499 //--------------------------------------------------------------------//
500 void AliPMDtracker::ResetClusters()
501 {
502   if (fRecpoints) fRecpoints->Clear();
503 }
504 //--------------------------------------------------------------------//