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