Correct and clean the vertex retrieval in case of SE or ME analysis
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliIsolationCut.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 /* $Id:  $ */
16
17 //_________________________________________________________________________
18 // Class containing methods for the isolation cut. 
19 // An AOD candidate (AliAODPWG4ParticleCorrelation type)
20 // is passed. Look in a cone around the candidate and study
21 // the hadronic activity inside to decide if the candidate is isolated
22 //
23 //
24 //*-- Author: Gustavo Conesa (LNF-INFN) 
25 //////////////////////////////////////////////////////////////////////////////
26   
27   
28 // --- ROOT system --- 
29 //#include <Riostream.h>
30 #include <TLorentzVector.h>
31 #include <TObjArray.h>
32
33 // --- AliRoot system --- 
34 #include "AliIsolationCut.h" 
35 #include "AliAODPWG4ParticleCorrelation.h"
36 #include "AliAODTrack.h"
37 #include "AliVCluster.h"
38 #include "AliCaloTrackReader.h"
39 #include "AliMixedEvent.h"
40
41 ClassImp(AliIsolationCut)
42   
43 //____________________________________________________________________________
44   AliIsolationCut::AliIsolationCut() : 
45     TObject(),
46     fConeSize(0.),fPtThreshold(0.), fPtFraction(0.), fICMethod(0),fPartInCone(0)
47  
48 {
49   //default ctor
50   
51   //Initialize parameters
52   InitParameters();
53
54 }
55 /*
56 //____________________________________________________________________________
57 AliIsolationCut::AliIsolationCut(const AliIsolationCut & g) : 
58   TObject(g),
59   fConeSize(g.fConeSize),
60   fPtThreshold(g.fPtThreshold),
61   fPtFraction(g.fPtFraction), 
62   fICMethod(g.fICMethod)
63 {
64   // cpy ctor
65   
66 }
67
68 //_________________________________________________________________________
69 AliIsolationCut & AliIsolationCut::operator = (const AliIsolationCut & source)
70 {
71   // assignment operator
72   
73   if(&source == this) return *this;
74   
75   fConeSize = source.fConeSize ;
76   fPtThreshold = source.fPtThreshold ; 
77   fICMethod = source.fICMethod ;
78   fPtFraction = source.fPtFraction ;
79
80   return *this;
81   
82 }
83 */
84 //____________________________________________________________________________
85 TString AliIsolationCut::GetICParametersList()
86 {
87   //Put data member values in string to keep in output container
88   
89   TString parList ; //this will be list of parameters used for this analysis.
90   const Int_t buffersize = 255;
91   char onePar[buffersize] ;
92   
93   snprintf(onePar,buffersize,"--- AliIsolationCut ---\n") ;
94   parList+=onePar ;     
95   snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
96   parList+=onePar ;
97   snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
98   parList+=onePar ;
99   snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
100   parList+=onePar ;
101   snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
102   parList+=onePar ;
103   snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
104   parList+=onePar ;
105
106   return parList; 
107 }
108
109 //____________________________________________________________________________
110 void AliIsolationCut::InitParameters()
111 {
112   //Initialize the parameters of the analysis.
113   
114   fConeSize    = 0.4 ; 
115   fPtThreshold = 1. ; 
116   fPtFraction  = 0.1 ; 
117   fPartInCone  = kNeutralAndCharged;
118   fICMethod    = kPtThresIC; // 0 pt threshol method, 1 cone pt sum method
119   
120 }
121
122 //__________________________________________________________________
123 void  AliIsolationCut::MakeIsolationCut(TObjArray * const plCTS,  TObjArray * const plNe, AliCaloTrackReader * const reader, 
124                                         const Bool_t fillAOD, AliAODPWG4ParticleCorrelation  *pCandidate, 
125                                         const TString aodArrayRefName,
126                                         Int_t & n, Int_t & nfrac, Float_t &coneptsum,  Bool_t  &isolated) const
127 {  
128   //Search in cone around a candidate particle if it is isolated 
129   Float_t phiC  = pCandidate->Phi() ;
130   Float_t etaC  = pCandidate->Eta() ;
131   Float_t ptC   = pCandidate->Pt() ;
132   Float_t pt    = -100. ;
133   Float_t eta   = -100.  ;
134   Float_t phi   = -100.  ;
135   Float_t rad   = -100 ;
136   n = 0 ;
137   coneptsum = 0.; 
138   isolated = kFALSE;
139   
140   //Initialize the array with refrences
141   TObjArray * refclusters = 0x0;
142   TObjArray * reftracks   = 0x0;
143   Int_t ntrackrefs   = 0;
144   Int_t nclusterrefs = 0;
145   
146   //Check charged particles in cone.
147   if(plCTS && (fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged)){
148     TVector3 p3;
149     for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
150       AliAODTrack* track = (AliAODTrack *)(plCTS->At(ipr)) ; 
151       //Do not count the candidate (pion, conversion photon) or the daughters of the candidate
152       if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1)) continue ;
153       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
154       pt   = p3.Pt();
155       eta  = p3.Eta();
156       phi  = p3.Phi() ;
157       if(phi<0) phi+=TMath::TwoPi();
158       
159       //Check if there is any particle inside cone with pt larger than  fPtThreshold
160       rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
161       
162       if(rad < fConeSize){
163         if(fillAOD) {
164           ntrackrefs++;
165           if(ntrackrefs == 1){
166             reftracks = new TObjArray(0);
167             reftracks->SetName(aodArrayRefName+"Tracks");
168             reftracks->SetOwner(kFALSE);
169           }
170           reftracks->Add(track);
171         }
172         //printf("charged in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
173         coneptsum+=pt;
174         if(pt > fPtThreshold ) n++;
175         if(pt > fPtFraction*ptC ) nfrac++;  
176       }
177     }// charged particle loop
178   }//Tracks
179   
180   //Check neutral particles in cone.  
181   if(plNe && (fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged)){
182           
183     //Get vertex for photon momentum calculation
184     //Double_t vertex2[] = {0,0,0} ; //vertex second AOD input ;
185     //if(reader->GetDataType()!= AliCaloTrackReader::kMC) 
186     //{
187       //if(reader->GetSecondInputAODTree()) reader->GetSecondInputAODVertex(vertex2);
188     //}
189     TLorentzVector mom ;
190     for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
191       AliVCluster * calo = (AliVCluster *)(plNe->At(ipr)) ;
192       
193       //Get the index where the cluster comes, to retrieve the corresponding vertex
194       Int_t evtIndex = 0 ; 
195       if (reader->GetMixedEvent()) {
196         evtIndex=reader->GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
197       }
198       
199       //Do not count the candidate (photon or pi0) or the daughters of the candidate
200       if(calo->GetID() == pCandidate->GetCaloLabel(0) || calo->GetID() == pCandidate->GetCaloLabel(1)) continue ;      //Skip matched clusters with tracks
201       
202       if(calo->GetNTracksMatched() > 0) continue ; 
203       
204       //Input from second AOD?
205       //Int_t input = 0;
206       //      if     (pCandidate->GetDetector() == "EMCAL" && reader->GetAODEMCALNormalInputEntries() <= ipr) input = 1 ;
207       //      else if(pCandidate->GetDetector() == "PHOS"  && reader->GetAODPHOSNormalInputEntries()  <= ipr) input = 1;
208       
209       //Get Momentum vector, 
210       //if     (input == 0) 
211       calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;//Assume that come from vertex in straight line
212       //else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line  
213       
214       pt   = mom.Pt();
215       eta  = mom.Eta();
216       phi  = mom.Phi() ;
217       if(phi<0) phi+=TMath::TwoPi();
218       
219       //Check if there is any particle inside cone with pt larger than  fPtThreshold
220       rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
221       if(rad < fConeSize){
222         if(fillAOD) {
223           nclusterrefs++;
224           if(nclusterrefs==1){
225             refclusters = new TObjArray(0);
226             refclusters->SetName(aodArrayRefName+"Clusters");
227             refclusters->SetOwner(kFALSE);
228           }
229           refclusters->Add(calo);
230         }
231         //printf("neutral in isolation cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
232         coneptsum+=pt;
233         if(pt > fPtThreshold ) n++;
234         if(pt > fPtFraction*ptC ) nfrac++;
235       }//in cone
236     }// neutral particle loop
237   }//neutrals
238   
239   //printf("Isolation Cut: in cone with: pT>pTthres %d, pT > pTfrac*pTcandidate %d \n",n,nfrac);
240   
241   //Add reference arrays to AOD when filling AODs only
242   if(fillAOD) {
243     if(refclusters)     pCandidate->AddObjArray(refclusters);
244     if(reftracks)         pCandidate->AddObjArray(reftracks);
245   }
246   
247   //Check isolation, depending on method.
248   if( fICMethod == kPtThresIC){
249     if(n==0) isolated = kTRUE ;
250   }
251   else if( fICMethod == kSumPtIC){
252     if(coneptsum < fPtThreshold)
253       isolated  =  kTRUE ;
254   }
255   else if( fICMethod == kPtFracIC){
256     if(nfrac==0) isolated = kTRUE ;
257   }
258   else if( fICMethod == kSumPtFracIC){
259     if(coneptsum < fPtFraction*ptC)
260       isolated  =  kTRUE ;
261   }
262   
263   //if(refclusters) delete refclusters;
264   //if(reftracks)   delete reftracks;
265   
266 }
267
268 //__________________________________________________________________
269 void AliIsolationCut::Print(const Option_t * opt) const
270 {
271   
272   //Print some relevant parameters set for the analysis
273   if(! opt)
274     return;
275   
276   printf("**** Print %s %s **** \n", GetName(), GetTitle() ) ;
277   
278   printf("IC method          =     %d\n", fICMethod) ; 
279   printf("Cone Size          =     %1.2f\n", fConeSize) ; 
280   printf("pT threshold       =     %2.1f\n", fPtThreshold) ;
281   printf("pT fraction        =     %3.1f\n", fPtFraction) ;
282   printf("particle type in cone =  %d\n",fPartInCone);
283   printf("    \n") ;
284   
285