TestBit for V0s
[u/mrichter/AliRoot.git] / PWGUD / UPC / AliAnalysisTaskUpcK0sK0s.cxx
CommitLineData
8b17ae4a
MB
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// c++ headers
17#include <iostream>
18#include <string.h>
86e51cbb 19#include <bitset>
8b17ae4a
MB
20
21// root headers
22#include "TH1I.h"
23#include "TTree.h"
24#include "TClonesArray.h"
25#include "TParticle.h"
18e1cc2e 26#include "TString.h"
8b17ae4a
MB
27#include "TFile.h"
28#include "TDatabasePDG.h"
29#include "TLorentzVector.h"
30
31// aliroot headers
32#include "AliAnalysisManager.h"
33#include "AliInputEventHandler.h"
34#include "AliESDEvent.h"
35#include "AliAODEvent.h"
36#include "AliMCEvent.h"
37#include "AliAODVZERO.h"
38#include "AliAODZDC.h"
39#include "AliESDVZERO.h"
40#include "AliESDZDC.h"
41#include "AliPIDResponse.h"
42#include "AliAODTrack.h"
43#include "AliAODPid.h"
44#include "AliAODVertex.h"
45#include "AliESDVertex.h"
46#include "AliMultiplicity.h"
47#include "AliESDtrack.h"
48#include "AliESDMuonTrack.h"
49#include "AliAODMCParticle.h"
50#include "AliMCParticle.h"
51
52// my headers
53#include "AliAnalysisTaskUpcK0sK0s.h"
54
55ClassImp(AliAnalysisTaskUpcK0sK0s);
56
57using std::cout;
58using std::endl;
59
60//trees for UPC analysis,
61// michal.broz@cern.ch
62
63//_____________________________________________________________________________
64AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s()
489951cd 65 : AliAnalysisTaskSE(),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
8b17ae4a
MB
66 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
67 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
68 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
69 fK0sAODv0s(0),fK0sESDv0s(0),fK0sAODTracks(0),fK0sESDTracks(0),
70 fListHist(0),fHistTriggersPerRun(0),fHistK0sMass(0)
71
72{
73
74//Dummy constructor
75
76}//AliAnalysisTaskUpcK0sK0s
77
78
79//_____________________________________________________________________________
80AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s(const char *name)
489951cd 81 : AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
8b17ae4a
MB
82 fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
83 fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
84 fDataFilnam(0),fRecoPass(0),fEvtNum(0),
85 fK0sAODv0s(0),fK0sESDv0s(0),fK0sAODTracks(0),fK0sESDTracks(0),
86 fListHist(0),fHistTriggersPerRun(0),fHistK0sMass(0)
87
88{
89
90 // Constructor
91 if( strstr(name,"ESD") ) fType = 0;
92 if( strstr(name,"AOD") ) fType = 1;
93
94 Init();
95
96 DefineOutput(1, TTree::Class());
489951cd 97 DefineOutput(2, TList::Class());
8b17ae4a
MB
98
99}//AliAnalysisTaskUpcK0sK0s
100
101//_____________________________________________________________________________
102void AliAnalysisTaskUpcK0sK0s::Init()
103{
104
105 for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
106
107}//Init
108
109//_____________________________________________________________________________
110AliAnalysisTaskUpcK0sK0s::~AliAnalysisTaskUpcK0sK0s()
111{
112 // Destructor
113 if(fK0sTree){
114 delete fK0sTree;
115 fK0sTree = 0x0;
116 }
8b17ae4a
MB
117 if(fListHist){
118 delete fListHist;
119 fListHist = 0x0;
120 }
121
122}//~AliAnalysisTaskUpcK0sK0s
123
124
125//_____________________________________________________________________________
126void AliAnalysisTaskUpcK0sK0s::UserCreateOutputObjects()
127{
8b17ae4a 128
8b17ae4a
MB
129 //vertices
130 fK0sAODv0s = new TClonesArray("AliAODv0", 1000);
131 fK0sESDv0s = new TClonesArray("AliESDv0", 1000);
132
133 //tracks
134 fK0sAODTracks = new TClonesArray("AliAODTrack", 1000);
135 fK0sESDTracks = new TClonesArray("AliESDtrack", 1000);
136
137 //output tree with K0s candidate events
138 fK0sTree = new TTree("fK0sTree", "fK0sTree");
139 fK0sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
140 fK0sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
141 fK0sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
142
143 fK0sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
144 fK0sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
145 fK0sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
146 fK0sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
147 fK0sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
148 fK0sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
149 fK0sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
150 fK0sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
151 fK0sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
152 fK0sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
153 fK0sTree ->Branch("fDataFilnam", &fDataFilnam);
154 fK0sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
155 fK0sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
156 if( fType == 0 ) {
157 fK0sTree ->Branch("fK0sESDv0s", &fK0sESDv0s);
158 fK0sTree ->Branch("fK0sESDTracks", &fK0sESDTracks);
159 }
160 if( fType == 1 ) {
161 fK0sTree ->Branch("fK0sAODv0s", &fK0sAODv0s);
162 fK0sTree ->Branch("fK0sAODTracks", &fK0sAODTracks);
163 }
164
165
166 fListHist = new TList();
167 fListHist ->SetOwner();
168
169 fHistTriggersPerRun = new TH1D("fHistTriggersPerRun", "fHistTriggersPerRun", 3000, 167000.5, 170000.5);
170 fListHist->Add(fHistTriggersPerRun);
171
489951cd 172 fHistK0sMass = new TH2D("fHistK0sMass","fHistK0sMass",20,0.45,0.55,20,0.45,0.55);
8b17ae4a
MB
173 fListHist->Add(fHistK0sMass);
174
175 PostData(1, fK0sTree);
489951cd 176 PostData(2, fListHist);
8b17ae4a
MB
177
178}//UserCreateOutputObjects
179
180//_____________________________________________________________________________
181void AliAnalysisTaskUpcK0sK0s::UserExec(Option_t *)
182{
183
184 //cout<<"#################### Next event ##################"<<endl;
185
186 if( fType == 0 ) RunESD();
187 if( fType == 1 ){
188 if(fRunHist) RunAODhist();
189 if(fRunTree) RunAODtree();
190 }
191
192}//UserExec
193//_____________________________________________________________________________
194void AliAnalysisTaskUpcK0sK0s::RunAODhist()
195{
18e1cc2e 196
8b17ae4a
MB
197 //input event
198 AliAODEvent *aod = (AliAODEvent*) InputEvent();
199 if(!aod) return;
86e51cbb 200
8b17ae4a
MB
201 //Trigger
202 TString trigger = aod->GetFiredTriggerClasses();
203
8a2486ba 204 if( !trigger.Contains("CCUP") ) return;
8b17ae4a
MB
205
206 fRunNum = aod ->GetRunNumber();
207 fHistTriggersPerRun->Fill(fRunNum);
208
8b17ae4a
MB
209 //primary vertex
210 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
211 fVtxContrib = fAODVertex->GetNContributors();
212 if(fVtxContrib < 2) return;
213
214
215 //VZERO, ZDC
216 AliAODVZERO *fV0data = aod ->GetVZEROData();
217
218 fV0Adecision = fV0data->GetV0ADecision();
219 fV0Cdecision = fV0data->GetV0CDecision();
220 if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
221
222 AliAODZDC *fZDCdata = aod->GetZDCData();
223
224 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
225 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
18e1cc2e 226 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
8b17ae4a
MB
227
228 Int_t nGoodV0s=0;
229 Int_t V0Index[3] = {-1,-1,-1};
230 Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
231
232 Int_t nGoodTracks=0;
233 Int_t TrackID[5] = {-1,-1,-1,-1,-1};
234
235
236 //V0s loop
237 for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
238 AliAODv0 *v0 = aod->GetV0(iV0);
239 if( !v0 ) continue;
240 Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
241 if (lOnFlyStatus) continue;
242
243 AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
244 AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
245 if (!pTrack || !nTrack) continue;
246
247 if ( pTrack->Charge() == nTrack->Charge())continue;
248
86e51cbb 249 if(!(pTrack->TestFilterBit(1<<0))) continue;
250 if(!(nTrack->TestFilterBit(1<<0))) continue;
18e1cc2e 251 if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
252 if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
8b17ae4a
MB
253 if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
254 if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
255 if(pTrack->GetTPCNcls() < 50)continue;
256 if(nTrack->GetTPCNcls() < 50)continue;
257 if(pTrack->Chi2perNDF() > 4)continue;
258 if(nTrack->Chi2perNDF() > 4)continue;
259
260 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
261 if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
262 if(TMath::Abs(dca[1]) > 2) continue;
263 if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
264 if(TMath::Abs(dca[1]) > 2) continue;
265
266 V0Index[nGoodV0s] = iV0;
267 V0TrackID[2*nGoodV0s] = pTrack->GetID();
268 V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
86e51cbb 269
8b17ae4a
MB
270 nGoodV0s++;
271 if(nGoodV0s > 2) break;
272 }//V0s loop
273
274 //Track loop
275 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
276 AliAODTrack *trk = aod->GetTrack(itr);
277 if( !trk ) continue;
8a2486ba 278 if(!(trk->TestFilterBit(1<<0))) continue;
8b17ae4a
MB
279
280 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
281 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
282 if(trk->GetTPCNcls() < 50)continue;
283 if(trk->Chi2perNDF() > 4)continue;
284 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
285 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
286 if(TMath::Abs(dca[1]) > 2) continue;
287
288 TrackID[nGoodTracks] = trk->GetID();
289 nGoodTracks++;
290
291 if(nGoodTracks > 4) break;
292 }//Track loop
293
294 if(nGoodV0s == 2 && nGoodTracks == 4){
86e51cbb 295 SortArray(TrackID);
296 SortArray(V0TrackID);
297 for(Int_t i=0; i<4; i++) if (TrackID[i] != V0TrackID[i]) return;
298
489951cd 299 AliAODv0 *v00 = aod->GetV0(V0Index[0]);
300 AliAODv0 *v01 = aod->GetV0(V0Index[1]);
18e1cc2e 301 fHistK0sMass->Fill(v00->MassK0Short(),v01->MassK0Short());
302
8b17ae4a 303 }
8b17ae4a 304
489951cd 305 PostData(2, fListHist);
8b17ae4a
MB
306
307}
308
309//_____________________________________________________________________________
310void AliAnalysisTaskUpcK0sK0s::RunAODtree()
311{
312 //input event
313 AliAODEvent *aod = (AliAODEvent*) InputEvent();
314 if(!aod) return;
315
316 //input data
18e1cc2e 317 fDataFilnam = aod->GetHeader()->GetESDFileName();
318 fEvtNum = aod->GetHeader()->GetEventNumberESDFile();
319 fRunNum = aod->GetRunNumber();
86e51cbb 320
8b17ae4a 321
8b17ae4a 322 //Trigger
86e51cbb 323 //Trigger
8b17ae4a
MB
324 TString trigger = aod->GetFiredTriggerClasses();
325
86e51cbb 326 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
327 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
328 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
329
330 Bool_t isTriggered = kFALSE;
331 for(Int_t i=0; i<ntrg; i++) {
332 if( fTrigger[i] ) isTriggered = kTRUE;
333 }
334 if(!isTriggered ) return;
8b17ae4a 335
8b17ae4a
MB
336 //trigger inputs
337 fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
338 fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
339
340 //Event identification
341 fPerNum = aod ->GetPeriodNumber();
342 fOrbNum = aod ->GetOrbitNumber();
343 fBCrossNum = aod ->GetBunchCrossNumber();
344
345 //primary vertex
346 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
347 fVtxContrib = fAODVertex->GetNContributors();
348
349 //Tracklets
350 fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
351
352 //VZERO, ZDC
353 AliAODVZERO *fV0data = aod ->GetVZEROData();
354 AliAODZDC *fZDCdata = aod->GetZDCData();
355
356 fV0Adecision = fV0data->GetV0ADecision();
357 fV0Cdecision = fV0data->GetV0CDecision();
358 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
359 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
360
361 Int_t nGoodV0s=0;
362 Int_t V0Index[3] = {-1,-1,-1};
363 Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
364
365 Int_t nGoodTracks=0;
366 Int_t TrackID[5] = {-1,-1,-1,-1,-1};
367
368 //V0s loop
369 for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
370 AliAODv0 *v0 = aod->GetV0(iV0);
371 if( !v0 ) continue;
372 Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
373 if (lOnFlyStatus) continue;
374
375 AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
376 AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
377 if (!pTrack || !nTrack) continue;
378
379 if ( pTrack->Charge() == nTrack->Charge())continue;
380
86e51cbb 381 if(!(pTrack->TestFilterBit(1<<0))) continue;
382 if(!(nTrack->TestFilterBit(1<<0))) continue;
8b17ae4a
MB
383 if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
384 if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
18e1cc2e 385 if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
386 if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
8b17ae4a
MB
387 if(pTrack->GetTPCNcls() < 50)continue;
388 if(nTrack->GetTPCNcls() < 50)continue;
389 if(pTrack->Chi2perNDF() > 4)continue;
390 if(nTrack->Chi2perNDF() > 4)continue;
391
392 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
393 if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
394 if(TMath::Abs(dca[1]) > 2) continue;
395 if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
396 if(TMath::Abs(dca[1]) > 2) continue;
397
398 V0Index[nGoodV0s] = iV0;
399 V0TrackID[2*nGoodV0s] = pTrack->GetID();
400 V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
401 nGoodV0s++;
402 if(nGoodV0s > 2) break;
403 }//V0s loop
404
405 //Track loop
406 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
407 AliAODTrack *trk = aod->GetTrack(itr);
408 if( !trk ) continue;
8a2486ba 409 if(!(trk->TestFilterBit(1<<0))) continue;
8b17ae4a
MB
410
411 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
412 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
413 if(trk->GetTPCNcls() < 50)continue;
414 if(trk->Chi2perNDF() > 4)continue;
415 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
416 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
417 if(TMath::Abs(dca[1]) > 2) continue;
418
419 TrackID[nGoodTracks] = trk->GetID();
420 nGoodTracks++;
421
422 if(nGoodTracks > 4) break;
423 }//Track loop
424
425 if(nGoodV0s == 2 && nGoodTracks == 4){
426 //SortArray(TrackID);
427 //SortArray(V0TrackID);
489951cd 428 //for{Int_t i=0; i<4; i++} if (TrackID[i] != V0TrackID[i]) return;
8b17ae4a
MB
429 for(Int_t i=0; i<2; i++){
430 AliAODv0 *v0 = aod->GetV0(V0Index[i]);
431 AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
432 AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
433
434 new((*fK0sAODv0s)[i]) AliAODv0(*v0);
435 new((*fK0sAODTracks)[2*i]) AliAODTrack(*pTrack);
436 new((*fK0sAODTracks)[2*i+1]) AliAODTrack(*nTrack);
437 }
438 fK0sTree ->Fill();
439
440 }
441
442 PostData(1, fK0sTree);
8b17ae4a
MB
443
444}//RunAOD
445
446//_____________________________________________________________________________
86e51cbb 447void AliAnalysisTaskUpcK0sK0s::SortArray(Int_t *dArray) {
8b17ae4a
MB
448 for (Int_t i = 3; i > 0; --i) {
449 for (Int_t j = 0; j < i; ++j) {
450 if (dArray[j] > dArray[j+1]) {
86e51cbb 451 Int_t dTemp = dArray[j];
8b17ae4a
MB
452 dArray[j] = dArray[j+1];
453 dArray[j+1] = dTemp;
454 }
455 }
456 }
457}
458
459//_____________________________________________________________________________
460void AliAnalysisTaskUpcK0sK0s::RunESD()
461{
462
8b17ae4a
MB
463}//RunESD
464
465//_____________________________________________________________________________
466void AliAnalysisTaskUpcK0sK0s::Terminate(Option_t *)
467{
468
469 cout<<"Analysis complete."<<endl;
470}//Terminate
471