Bug fix , problem was with pile-up tagging
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / MultEvShape / AliAnalysisCentralCutEvtESD.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2008, 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 //  ESD event level cuts for azimuthal isotropic
18 //  expansion in highly central collisions analysis
19 //  author: Cristian Andrei
20 //           acristian@niham.nipne.ro
21 //  --------------------------------------------------
22
23 #include "TMath.h"
24 #include <TObjArray.h>
25 #include <TList.h>
26
27 #include "AliESDEvent.h"
28 #include "AliESDtrack.h"
29 #include "AliESDtrackCuts.h"
30  #include "AliMultiplicity.h"
31
32 #include "AliAnalysisCentralCutESD.h"
33 #include "AliAnalysisCentralCutEvtESD.h"
34
35 class TObject;
36
37
38 //____________________________________________________________________
39 ClassImp(AliAnalysisCentralCutEvtESD)
40
41 //____________________________________________________________________
42 AliAnalysisCentralCutEvtESD::AliAnalysisCentralCutEvtESD(const Char_t* name, const Char_t* title) 
43     :AliAnalysisCuts(name,title)
44     ,fReqMult(kFALSE)
45     ,fReqDir(kFALSE)
46     ,fReqSPDMult(kFALSE)
47     ,fReqSPDDir(kFALSE)
48     ,fMultMin(0)
49     ,fMultMax(0)
50     ,fDirMin(0)
51     ,fDirMax(0)
52     ,fSPDMultMin(0)
53     ,fSPDMultMax(0)
54     ,fSPDDirMin(0)
55     ,fSPDDirMax(0)
56
57 {
58 //constructor
59     for(Int_t i=0; i<10; i++){
60                 fCutsList[i] = 0;
61     }
62
63         InitCuts();
64
65 }
66
67 AliAnalysisCentralCutEvtESD::~AliAnalysisCentralCutEvtESD() {
68 //Destructor
69   for (Int_t i=0; i<10; i++)
70     delete fCutsList[i];
71
72 }
73
74
75
76 Bool_t AliAnalysisCentralCutEvtESD::IsSelected(TObject *obj){
77 // check whether the event passes the cuts
78     AliESDEvent *esdEvent = dynamic_cast<AliESDEvent *>(obj);
79     if(!esdEvent){
80                 printf("AliAnalysisCentralCutEvtESD:IsSelected ->Can't get ESD event!\n");
81                 return kFALSE;
82     }
83
84         if(fReqMult){
85                 Int_t mult = CalcMult(esdEvent);
86                 if((mult<fMultMin)||(mult>fMultMax)){
87                         return kFALSE;
88                 }
89     }
90
91     if(fReqDir){
92                 Double_t dir = CalcDir(esdEvent);
93                 if((dir<fDirMin)||(dir>fDirMax)){
94                         return kFALSE;
95                 }
96     }
97
98     if(fReqSPDMult){
99                 Double_t spdMult = CalcSPDMult(esdEvent);
100                 if((spdMult<fSPDMultMin)||(spdMult>fSPDMultMax)){
101                         return kFALSE;
102                 }
103     }
104
105     if(fReqSPDDir){
106                 Double_t spdDir = CalcSPDDir(esdEvent);
107                 if((spdDir<fSPDDirMin)||(spdDir>fSPDDirMax)){
108                         return kFALSE;
109                 }
110     }
111
112     return kTRUE;
113 }
114
115 //_________________________________________________________________________
116 void AliAnalysisCentralCutEvtESD::InitCuts(){
117 //Initialize internal cuts
118
119 ////////////////  FULL ALICE ///////////////////
120 //------------General ESD Cuts------------------
121     AliESDtrackCuts *esdCutsGen = new AliESDtrackCuts("AliESDtrackCuts", "CutEvtESDInternal");
122     esdCutsGen->SetMinNClustersTPC(50);
123     esdCutsGen->SetMaxChi2PerClusterTPC(2.2);
124     esdCutsGen->SetMaxCovDiagonalElements(0.5,0.5,0.5,0.5,0.5);
125     esdCutsGen->SetRequireTPCRefit(kTRUE);
126     esdCutsGen->SetAcceptKinkDaughters(kFALSE);
127     esdCutsGen->SetMaxNsigmaToVertex(2.0);
128     esdCutsGen->SetRequireSigmaToVertex(kTRUE);
129
130     AliAnalysisCentralCutESD *esdCutsGen1 = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
131     esdCutsGen1->SetReqIsCharged();
132
133 //-------------Specific ESD Cuts------------------
134     AliESDtrackCuts *esdCutsMult = new AliESDtrackCuts("AliESDCutsMult", "midrapidity");
135     esdCutsMult->SetEtaRange(-0.5,0.5);
136     AliESDtrackCuts *esdCutsDir = new AliESDtrackCuts("AliESDCutsDir", "SPD coverage");
137     esdCutsDir->SetEtaRange(0.0,1.9);
138
139 ///////////////////  SPD ONLY  ////////////////////
140     AliESDtrackCuts *esdCutsSPD = new AliESDtrackCuts("AliESDtrackCuts", "CutEvtESDInternal");
141     esdCutsSPD->SetAcceptKinkDaughters(kFALSE);
142     esdCutsSPD->SetMaxNsigmaToVertex(2.0);
143     esdCutsSPD->SetRequireSigmaToVertex(kTRUE);
144
145 //--------------set the cuts ----------------------
146     TObjArray* esdListMult = new TObjArray();
147     esdListMult->AddLast(esdCutsGen);
148     esdListMult->AddLast(esdCutsGen1);
149     esdListMult->AddLast(esdCutsMult);
150
151     TObjArray* esdListDir = new TObjArray();
152     esdListDir->AddLast(esdCutsGen);
153     esdListDir->AddLast(esdCutsGen1);
154     esdListDir->AddLast(esdCutsDir);
155
156     TObjArray* esdListSPD = new TObjArray();
157     esdListSPD->AddLast(esdCutsSPD);
158
159     fCutsList[0]=esdListDir;
160     fCutsList[1]=esdListMult;
161     fCutsList[2]=esdListSPD;
162
163 }
164
165
166 //__________________________________________________________________________
167 Bool_t AliAnalysisCentralCutEvtESD::CheckIntCuts(Int_t no, TObject *obj) const{
168 // Check if the particle passes the internal cuts
169
170     if(no > 9){
171                 printf("\n AliAnalysisCentralCutEvtESD::CheckIntCuts -> Cut number is not ok! \n");
172                 return kFALSE;
173     }
174
175     if(!fCutsList[no]){
176                 printf("AliAnalysisCentralCutEvtESD::CheckIntCuts -> cuts list problem! \n");
177                 return kFALSE;
178     }
179
180
181     TObjArrayIter iter(fCutsList[no]);
182     AliAnalysisCuts *cut = 0;
183
184
185     while((cut = (AliAnalysisCuts*)iter.Next())){
186
187                 if(!cut->IsSelected(obj)) return kFALSE;
188     }
189
190     return kTRUE;
191 }
192
193
194 //__________________________________________________________________________
195 Double_t AliAnalysisCentralCutEvtESD::CalcDir(AliESDEvent* const esdEv) {
196
197 //Compute the directivity - FULL ALICE
198
199     Double_t dir;
200     Double_t px,py;
201     Double_t sumaPt = 0;
202     Double_t sumaPx = 0;
203     Double_t sumaPy = 0;
204
205     Double_t pt;
206
207     if (!esdEv){
208                 printf("NULL tree\n");
209                 return -1;
210     }
211
212     Int_t totTracks=esdEv->GetNumberOfTracks();
213
214     for(Int_t itrack = 0; itrack < totTracks; itrack++){//first loop->compute event directivity
215
216                 AliESDtrack* track = esdEv->GetTrack(itrack);
217                 
218                 if (!track) {
219                                 Printf("ERROR: Could not receive track %d", itrack);
220                         continue;
221                 }
222         
223                 if(!CheckIntCuts(0, track)) continue;
224                 
225                 px = track->Px();       
226                 py = track->Py();       
227                 pt = track->Pt();
228         
229                 sumaPx = sumaPx + px;
230                 sumaPy = sumaPy + py;
231         
232                 sumaPt = sumaPt + pt;
233
234     }//end track loop
235
236
237         if(sumaPt < 0.0000001){
238                 return -1;
239     }
240
241     dir = (sqrt(pow(sumaPx,2)+pow(sumaPy,2)))/sumaPt;
242
243     return dir;
244 }
245
246 //__________________________________________________________________________
247 Int_t AliAnalysisCentralCutEvtESD::CalcMult(AliESDEvent* const esdEv) {
248
249 //Compute multiplicity - FULL ALICE
250
251     Int_t charged = 0;
252
253     if (!esdEv){
254                 printf("NULL tree\n");
255                 return -1;
256     }
257
258     Int_t totTracks=esdEv->GetNumberOfTracks();
259
260     for(Int_t iTrack = 0; iTrack < totTracks; iTrack++){//second track loop -> compute event multiplicity
261
262                 AliESDtrack* track = esdEv->GetTrack(iTrack);
263                 
264                 if (!track) {
265                         Printf("ERROR: Could not receive track %d", iTrack);
266                         continue;
267                 }
268         
269                 if(!CheckIntCuts(1, track)) continue;
270         
271                 
272                 charged++; //multiplicity
273
274     }//end second track loop
275
276
277     return charged;
278 }
279
280
281 //__________________________________________________________________________
282 Double_t AliAnalysisCentralCutEvtESD::CalcSPDDir(AliESDEvent* const esdEv) {
283
284 //Compute directivity - SPD ONLY
285
286     Double_t dirU;
287     Double_t pxU,pyU;
288     Double_t sumaPxU = 0;
289     Double_t sumaPyU = 0;
290
291     Double_t goodtrack = 0;
292     Double_t spdEta = 0.;
293
294     const AliMultiplicity *spdMult=esdEv->GetMultiplicity();
295     if(!spdMult){
296         printf("Unable to get multiplicity! \n");
297         return -1;
298     }
299
300     Int_t spdTracks = spdMult->GetNumberOfTracklets();  //SPD multiplicity
301
302     for(Int_t iTrack = 0; iTrack< spdTracks; iTrack++){  //SPD track loop -> Directivity
303         
304                 AliESDtrack* track = esdEv->GetTrack(iTrack);
305                 
306                 if (!track) {
307                                 Printf("ERROR: Could not receive track %d", iTrack);
308                                 continue;
309                 }
310         
311                 if(!CheckIntCuts(2, track)) continue; 
312                 
313                 spdEta = spdMult->GetEta(iTrack);
314                 
315                 if((spdEta<0.0)||(spdEta>1.9)) continue;
316         
317                 Double_t phi = spdMult->GetPhi(iTrack);  
318         
319                 pxU = TMath::Cos(phi);
320                 pyU = TMath::Sin(phi);
321         
322                 sumaPxU = sumaPxU + pxU;
323                 sumaPyU = sumaPyU + pyU;        
324                 
325                 goodtrack++;
326         
327     }//end SPD track loop
328
329     if(goodtrack < 1.) return -1;
330
331     dirU = (sqrt(pow(sumaPxU,2)+pow(sumaPyU,2)))/goodtrack;
332
333     return dirU;
334 }
335
336 //__________________________________________________________________________
337
338 Int_t AliAnalysisCentralCutEvtESD::CalcSPDMult(AliESDEvent* const esdEv) {
339
340         //Compute multiplicity - SPD ONLY
341
342     const AliMultiplicity *spdMult=esdEv->GetMultiplicity();
343     if(!spdMult){
344         printf("Unable to get multiplicity! \n");
345         return -1;
346     }
347
348     Int_t spdTracks = spdMult->GetNumberOfTracklets();  //SPD multiplicity
349
350     Double_t spdEta;
351     Int_t charged = 0;
352
353     for(Int_t iTrack = 0; iTrack< spdTracks; iTrack++){  //second SPD track loop -> Multiplicity
354         
355                 AliESDtrack* track = esdEv->GetTrack(iTrack);
356                 
357                 if (!track) {
358                                 Printf("ERROR: Could not receive track %d", iTrack);
359                                 continue;
360                 }
361                 
362                 if(!CheckIntCuts(2, track)) continue;
363                 
364                 spdEta = spdMult->GetEta(iTrack);
365                 
366                 if((spdEta<-0.5)||(spdEta>0.5)) continue;
367         
368                 charged++;
369         
370     }//end second SPD track loop 
371
372     return charged;
373 }