1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
24 #include "TClonesArray.h"
25 #include "TParticle.h"
28 #include "TDatabasePDG.h"
29 #include "TLorentzVector.h"
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"
53 #include "AliAnalysisTaskUpcK0sK0s.h"
55 ClassImp(AliAnalysisTaskUpcK0sK0s);
60 //trees for UPC analysis,
61 // michal.broz@cern.ch
63 //_____________________________________________________________________________
64 AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s()
65 : AliAnalysisTaskSE(),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
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)
76 }//AliAnalysisTaskUpcK0sK0s
79 //_____________________________________________________________________________
80 AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s(const char *name)
81 : AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
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)
91 if( strstr(name,"ESD") ) fType = 0;
92 if( strstr(name,"AOD") ) fType = 1;
96 DefineOutput(1, TTree::Class());
97 DefineOutput(2, TList::Class());
99 }//AliAnalysisTaskUpcK0sK0s
101 //_____________________________________________________________________________
102 void AliAnalysisTaskUpcK0sK0s::Init()
105 for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
109 //_____________________________________________________________________________
110 AliAnalysisTaskUpcK0sK0s::~AliAnalysisTaskUpcK0sK0s()
122 }//~AliAnalysisTaskUpcK0sK0s
125 //_____________________________________________________________________________
126 void AliAnalysisTaskUpcK0sK0s::UserCreateOutputObjects()
130 fK0sAODv0s = new TClonesArray("AliAODv0", 1000);
131 fK0sESDv0s = new TClonesArray("AliESDv0", 1000);
134 fK0sAODTracks = new TClonesArray("AliAODTrack", 1000);
135 fK0sESDTracks = new TClonesArray("AliESDtrack", 1000);
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");
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");
157 fK0sTree ->Branch("fK0sESDv0s", &fK0sESDv0s);
158 fK0sTree ->Branch("fK0sESDTracks", &fK0sESDTracks);
161 fK0sTree ->Branch("fK0sAODv0s", &fK0sAODv0s);
162 fK0sTree ->Branch("fK0sAODTracks", &fK0sAODTracks);
166 fListHist = new TList();
167 fListHist ->SetOwner();
169 fHistTriggersPerRun = new TH1D("fHistTriggersPerRun", "fHistTriggersPerRun", 3000, 167000.5, 170000.5);
170 fListHist->Add(fHistTriggersPerRun);
172 fHistK0sMass = new TH2D("fHistK0sMass","fHistK0sMass",20,0.45,0.55,20,0.45,0.55);
173 fListHist->Add(fHistK0sMass);
175 PostData(1, fK0sTree);
176 PostData(2, fListHist);
178 }//UserCreateOutputObjects
180 //_____________________________________________________________________________
181 void AliAnalysisTaskUpcK0sK0s::UserExec(Option_t *)
184 //cout<<"#################### Next event ##################"<<endl;
186 if( fType == 0 ) RunESD();
188 if(fRunHist) RunAODhist();
189 if(fRunTree) RunAODtree();
193 //_____________________________________________________________________________
194 void AliAnalysisTaskUpcK0sK0s::RunAODhist()
198 AliAODEvent *aod = (AliAODEvent*) InputEvent();
202 TString trigger = aod->GetFiredTriggerClasses();
204 if( !trigger.Contains("CCUP") ) return;
206 fRunNum = aod ->GetRunNumber();
207 fHistTriggersPerRun->Fill(fRunNum);
210 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
211 fVtxContrib = fAODVertex->GetNContributors();
212 if(fVtxContrib < 2) return;
216 AliAODVZERO *fV0data = aod ->GetVZEROData();
218 fV0Adecision = fV0data->GetV0ADecision();
219 fV0Cdecision = fV0data->GetV0CDecision();
220 if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
222 AliAODZDC *fZDCdata = aod->GetZDCData();
224 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
225 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
226 if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
229 Int_t V0Index[3] = {-1,-1,-1};
230 Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
233 Int_t TrackID[5] = {-1,-1,-1,-1,-1};
237 for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
238 AliAODv0 *v0 = aod->GetV0(iV0);
240 Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
241 if (lOnFlyStatus) continue;
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;
247 if ( pTrack->Charge() == nTrack->Charge())continue;
249 if(!(pTrack->TestFilterBit(1<<0))) continue;
250 if(!(nTrack->TestFilterBit(1<<0))) continue;
251 if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
252 if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
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;
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;
266 V0Index[nGoodV0s] = iV0;
267 V0TrackID[2*nGoodV0s] = pTrack->GetID();
268 V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
271 if(nGoodV0s > 2) break;
275 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
276 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
277 if(!trk) AliFatal("Not a standard AOD");
279 if(!(trk->TestFilterBit(1<<0))) continue;
281 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
282 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
283 if(trk->GetTPCNcls() < 50)continue;
284 if(trk->Chi2perNDF() > 4)continue;
285 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
286 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
287 if(TMath::Abs(dca[1]) > 2) continue;
289 TrackID[nGoodTracks] = trk->GetID();
292 if(nGoodTracks > 4) break;
295 if(nGoodV0s == 2 && nGoodTracks == 4){
297 SortArray(V0TrackID);
298 for(Int_t i=0; i<4; i++) if (TrackID[i] != V0TrackID[i]) return;
300 AliAODv0 *v00 = aod->GetV0(V0Index[0]);
301 AliAODv0 *v01 = aod->GetV0(V0Index[1]);
302 fHistK0sMass->Fill(v00->MassK0Short(),v01->MassK0Short());
306 PostData(2, fListHist);
310 //_____________________________________________________________________________
311 void AliAnalysisTaskUpcK0sK0s::RunAODtree()
314 AliAODEvent *aod = (AliAODEvent*) InputEvent();
318 AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
319 if(!header) AliFatal("Not a standard AOD");
321 fDataFilnam = header->GetESDFileName();
322 fEvtNum = header->GetEventNumberESDFile();
323 fRunNum = aod->GetRunNumber();
328 TString trigger = aod->GetFiredTriggerClasses();
330 fTrigger[0] = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
331 fTrigger[1] = trigger.Contains("CCUP2-B"); // Double gap
332 fTrigger[2] = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
334 Bool_t isTriggered = kFALSE;
335 for(Int_t i=0; i<ntrg; i++) {
336 if( fTrigger[i] ) isTriggered = kTRUE;
338 if(!isTriggered ) return;
341 fL0inputs = ((AliVAODHeader*)aod->GetHeader())->GetL0TriggerInputs();
342 fL1inputs = ((AliVAODHeader*)aod->GetHeader())->GetL1TriggerInputs();
344 //Event identification
345 fPerNum = aod ->GetPeriodNumber();
346 fOrbNum = aod ->GetOrbitNumber();
347 fBCrossNum = aod ->GetBunchCrossNumber();
350 AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
351 fVtxContrib = fAODVertex->GetNContributors();
354 fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
357 AliAODVZERO *fV0data = aod ->GetVZEROData();
358 AliAODZDC *fZDCdata = aod->GetZDCData();
360 fV0Adecision = fV0data->GetV0ADecision();
361 fV0Cdecision = fV0data->GetV0CDecision();
362 fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
363 fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
366 Int_t V0Index[3] = {-1,-1,-1};
367 Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
370 Int_t TrackID[5] = {-1,-1,-1,-1,-1};
373 for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
374 AliAODv0 *v0 = aod->GetV0(iV0);
376 Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
377 if (lOnFlyStatus) continue;
379 AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
380 AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
381 if (!pTrack || !nTrack) continue;
383 if ( pTrack->Charge() == nTrack->Charge())continue;
385 if(!(pTrack->TestFilterBit(1<<0))) continue;
386 if(!(nTrack->TestFilterBit(1<<0))) continue;
387 if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
388 if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
389 if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
390 if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
391 if(pTrack->GetTPCNcls() < 50)continue;
392 if(nTrack->GetTPCNcls() < 50)continue;
393 if(pTrack->Chi2perNDF() > 4)continue;
394 if(nTrack->Chi2perNDF() > 4)continue;
396 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
397 if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
398 if(TMath::Abs(dca[1]) > 2) continue;
399 if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
400 if(TMath::Abs(dca[1]) > 2) continue;
402 V0Index[nGoodV0s] = iV0;
403 V0TrackID[2*nGoodV0s] = pTrack->GetID();
404 V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
406 if(nGoodV0s > 2) break;
410 for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
411 AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
412 if(!trk) AliFatal("Not a standard AOD");
414 if(!(trk->TestFilterBit(1<<0))) continue;
416 if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
417 if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
418 if(trk->GetTPCNcls() < 50)continue;
419 if(trk->Chi2perNDF() > 4)continue;
420 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
421 if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
422 if(TMath::Abs(dca[1]) > 2) continue;
424 TrackID[nGoodTracks] = trk->GetID();
427 if(nGoodTracks > 4) break;
430 if(nGoodV0s == 2 && nGoodTracks == 4){
431 //SortArray(TrackID);
432 //SortArray(V0TrackID);
433 //for{Int_t i=0; i<4; i++} if (TrackID[i] != V0TrackID[i]) return;
434 for(Int_t i=0; i<2; i++){
435 AliAODv0 *v0 = aod->GetV0(V0Index[i]);
436 AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
437 AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
439 new((*fK0sAODv0s)[i]) AliAODv0(*v0);
440 new((*fK0sAODTracks)[2*i]) AliAODTrack(*pTrack);
441 new((*fK0sAODTracks)[2*i+1]) AliAODTrack(*nTrack);
447 PostData(1, fK0sTree);
451 //_____________________________________________________________________________
452 void AliAnalysisTaskUpcK0sK0s::SortArray(Int_t *dArray) {
453 for (Int_t i = 3; i > 0; --i) {
454 for (Int_t j = 0; j < i; ++j) {
455 if (dArray[j] > dArray[j+1]) {
456 Int_t dTemp = dArray[j];
457 dArray[j] = dArray[j+1];
464 //_____________________________________________________________________________
465 void AliAnalysisTaskUpcK0sK0s::RunESD()
470 //_____________________________________________________________________________
471 void AliAnalysisTaskUpcK0sK0s::Terminate(Option_t *)
474 cout<<"Analysis complete."<<endl;