Updating track cuts and removing functionality that is now in separate cosmics tasks
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliPWG4CosmicCandidates.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 // Analysis task looking for cosmic candidates embedded in a p-p event
18 // Task checks if particles are back-to-back in eta and phi
19 //             
20 //
21 // Author : Marta Verweij - UU - marta.verweij@cern.ch
22 //-----------------------------------------------------------------------
23
24 #include "TVector3.h"
25 #include <iostream>
26 #include "TH1.h"
27 #include "TH2.h"
28 #include "TH3.h"
29 #include "TList.h"
30 #include "TChain.h"
31
32 #include "AliAnalysisTask.h"
33 #include "AliAnalysisManager.h"
34
35 #include "AliESDEvent.h"
36 #include "AliESDInputHandler.h"
37
38 #include "AliESDtrack.h"
39 #include "AliESDtrackCuts.h"
40 #include "AliExternalTrackParam.h"
41
42 #include "AliLog.h"
43
44 #include "AliPWG4CosmicCandidates.h"
45
46 //using namespace std; //required for resolving the 'cout' symbol
47 using namespace std;
48
49 ClassImp(AliPWG4CosmicCandidates)
50
51 //________________________________________________________________________
52 AliPWG4CosmicCandidates::AliPWG4CosmicCandidates()
53 : AliAnalysisTaskSE(),
54   fTrackCuts(0), 
55   fPtMin(5.),
56   fMaxCosmicAngle(0.002),
57   fNEventAll(0),
58   fNEventSel(0),
59   fPtSignedCosmicCandidates(0),
60   fDeltaPtCosmicCandidates(0),
61   fDeltaPhiSumEta(0),
62   fDCAZCosmicCandidates(0),
63   fDCARCosmicCandidates(0),
64   fTheta(0),
65   fThetaZoom(0),
66   fThetaPt1Pt2(0),
67   fThetaPt1Pt2Signed(0),
68   fDeltaPhiSumEtaPt1(0),
69   fDeltaPhiSumEtaPt2(0),
70   fThetaDCAZ1DCAZ2(0),
71   fRisol(0),
72   fRisolTheta(0),
73   fHistListCosmics(0)
74 {
75   //
76   // Default constructor
77   //
78 }
79
80 //________________________________________________________________________
81 AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const char *name)
82   : AliAnalysisTaskSE(name),
83     fTrackCuts(0), 
84     fPtMin(5.),
85     fMaxCosmicAngle(0.002),
86     fNEventAll(0),
87     fNEventSel(0),
88     fPtSignedCosmicCandidates(0),
89     fDeltaPtCosmicCandidates(0),
90     fDeltaPhiSumEta(0),
91     fDCAZCosmicCandidates(0),
92     fDCARCosmicCandidates(0),
93     fTheta(0),
94     fThetaZoom(0),
95     fThetaPt1Pt2(0),
96     fThetaPt1Pt2Signed(0),
97     fDeltaPhiSumEtaPt1(0),
98     fDeltaPhiSumEtaPt2(0),
99     fThetaDCAZ1DCAZ2(0),
100     fRisol(0),
101     fRisolTheta(0),
102     fHistListCosmics(0)
103 {
104   // Constructor. Initialization of Inputs and Outputs
105
106   // Define input and output slots here
107   // Input slot #0 works with a TChain
108   //  DefineInput(0, TChain::Class());
109   // Output slot #0 id reserved by the base class for AOD
110   // Output slot #1 writes into a TList
111   DefineOutput(1, TList::Class());
112   // Output slot #2 writes into a AliESDtrackCuts
113   DefineOutput(2, AliESDtrackCuts::Class());
114 }
115 //________________________________________________________________________
116 AliPWG4CosmicCandidates::AliPWG4CosmicCandidates(const AliPWG4CosmicCandidates &res)
117   : AliAnalysisTaskSE(res),
118     fTrackCuts(0), 
119     fPtMin(5.),
120     fMaxCosmicAngle(0.002),
121     fNEventAll(0),
122     fNEventSel(0),
123     fPtSignedCosmicCandidates(0),
124     fDeltaPtCosmicCandidates(0),
125     fDeltaPhiSumEta(0),
126     fDCAZCosmicCandidates(0),
127     fDCARCosmicCandidates(0),
128     fTheta(0),
129     fThetaZoom(0),
130     fThetaPt1Pt2(0),
131     fThetaPt1Pt2Signed(0),
132     fDeltaPhiSumEtaPt1(0),
133     fDeltaPhiSumEtaPt2(0),
134     fThetaDCAZ1DCAZ2(0),
135     fRisol(0),
136     fRisolTheta(0),
137     fHistListCosmics(0)
138 {
139   //
140   // Dummy copy constructor
141   //
142 }
143
144 //________________________________________________________________________
145 AliPWG4CosmicCandidates& AliPWG4CosmicCandidates::operator=(const AliPWG4CosmicCandidates& /*trclass*/)
146 {
147   //
148   // Dummy assignment operator
149   //
150   return *this;
151 }
152
153 //________________________________________________________________________
154 void AliPWG4CosmicCandidates::LocalInit()
155 {
156   //
157   // Only called once at beginning
158   //
159   PostData(2,fTrackCuts);
160 }
161
162 //________________________________________________________________________
163 void AliPWG4CosmicCandidates::UserCreateOutputObjects()
164 {
165   //
166   // Create output objects
167   // Called once
168   AliDebug(2,Form(">> AliPWG4CosmicCandidates::UserCreateOutputObjects \n")); 
169
170   Bool_t oldStatus = TH1::AddDirectoryStatus();
171   TH1::AddDirectory(kFALSE); 
172   
173   OpenFile(1);
174   fHistListCosmics = new TList();
175
176   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
177   fHistListCosmics->Add(fNEventAll);
178   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
179   fHistListCosmics->Add(fNEventSel);
180
181   Float_t fgkPtMin=0.;
182   Float_t fgkPtMax=100.;
183   Int_t fgkNPtBins= (int)(fgkPtMax-fgkPtMin);
184
185   Int_t fgkNPhiBins=18;
186   Float_t kMinPhi = -0.5*TMath::Pi();
187   Float_t kMaxPhi = 3./2.*TMath::Pi();
188
189   Int_t fgkNThetaBins=fgkNPhiBins*8;
190   Float_t kMinTheta = -0.5*TMath::Pi();
191   Float_t kMaxTheta = 3./2.*TMath::Pi();
192
193   Int_t fgkNDCARBins=80;
194   Float_t fgkDCARMin = -0.2;
195   Float_t fgkDCARMax = 0.2;
196   Float_t *binsDCAR=new Float_t[fgkNDCARBins+1];
197   for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Float_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Float_t)i ;
198
199   Int_t fgkNDCAZBins=80;
200   Float_t fgkDCAZMin = -2.;
201   Float_t fgkDCAZMax = 2.;
202   Float_t *binsDCAZ=new Float_t[fgkNDCAZBins+1];
203   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Float_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Float_t)i ;
204
205   fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2*(int)(fgkPtMax-fgkPtMin), -1.*fgkPtMax, fgkPtMax);
206   fHistListCosmics->Add(fPtSignedCosmicCandidates);  
207
208   fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
209   fHistListCosmics->Add(fDeltaPtCosmicCandidates);  
210
211   fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,kMinPhi,kMaxPhi,80, -2.,2.);
212   fHistListCosmics->Add(fDeltaPhiSumEta);  
213
214   fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
215   fHistListCosmics->Add(fDCAZCosmicCandidates);
216
217   fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
218   fHistListCosmics->Add(fDCARCosmicCandidates);
219
220   fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
221   fHistListCosmics->Add(fTheta);
222
223   fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
224   fHistListCosmics->Add(fThetaZoom);
225
226   fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
227   fHistListCosmics->Add(fThetaPt1Pt2);
228
229   fThetaPt1Pt2Signed = new TH3F("fThetaPt1Pt2Signed","fThetaPt1Pt2Signed",fgkNThetaBins,kMinTheta,kMaxTheta,4*(int)(fgkPtMax-fgkPtMin),-1.*fgkPtMax,fgkPtMax,4*(int)(fgkPtMax-fgkPtMin),-1.*fgkPtMax,fgkPtMax);
230   fHistListCosmics->Add(fThetaPt1Pt2Signed);
231
232   fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
233   fHistListCosmics->Add(fDeltaPhiSumEtaPt1);
234
235   fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
236   fHistListCosmics->Add(fDeltaPhiSumEtaPt2);
237
238   fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
239   fHistListCosmics->Add(fThetaDCAZ1DCAZ2);
240
241   fRisol = new TH1F("fRisol","fRisol",100,0.,10.);
242   fHistListCosmics->Add(fRisol);
243
244   fRisolTheta = new TH2F("fRisolTheta","fRisolTheta",100,0.,10.,fgkNThetaBins,kMinTheta,kMaxTheta);
245   fHistListCosmics->Add(fRisolTheta);
246
247   TH1::AddDirectory(oldStatus);   
248
249   if(binsDCAR) delete [] binsDCAR;
250   if(binsDCAZ) delete [] binsDCAZ;
251
252 }
253
254 //________________________________________________________________________
255 void AliPWG4CosmicCandidates::UserExec(Option_t *) 
256 {
257   // Main loop
258   // Called for each event
259   //
260
261   // All events without selection
262   fNEventAll->Fill(0.);
263    
264   if (!fInputEvent) {
265     AliDebug(2,Form("ERROR: fESD not available"));
266     cout << "ERROR: fESD not available" << endl;
267     PostData(1, fHistListCosmics);
268     return;
269   }
270
271   const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
272   // Need vertex cut
273   if (!vtx || vtx->GetNContributors() < 2) {
274     // Post output data
275     PostData(1, fHistListCosmics);
276     return;
277   }
278
279   //  AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
280   double primVtx[3];
281   vtx->GetXYZ(primVtx);
282   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
283     // Post output data
284     PostData(1, fHistListCosmics);
285     return;
286   }
287   if(!fInputEvent->GetNumberOfTracks() || fInputEvent->GetNumberOfTracks()<2){ 
288     // Post output data
289     PostData(1, fHistListCosmics);
290     return;
291   }
292   Int_t nTracks = fInputEvent->GetNumberOfTracks();
293
294   if(!fTrackCuts) {
295    // Post output data
296     PostData(1, fHistListCosmics);
297     return;
298   }
299
300   fNEventSel->Fill(0.);
301
302   Float_t dcaR[2] = {0.,0.};
303   Float_t dcaZ[2] = {0.,0.};
304
305   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
306
307     AliESDtrack* track1 = (AliESDtrack*)fInputEvent->GetTrack(iTrack1);
308     if (!track1)  continue;
309     if(!(fTrackCuts->AcceptTrack(track1))) { continue; }
310     if(track1->Pt()<fPtMin) continue;
311     //Start 2nd track loop to look for correlations
312     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
313       AliESDtrack *track2 = (AliESDtrack*)fInputEvent->GetTrack(iTrack2);
314       if(!track2) continue;
315       if(!(fTrackCuts->AcceptTrack(track2))) { continue; }
316          
317       //Check if back-to-back
318       Double_t mom1[3],mom2[3];
319       track1->GetPxPyPz(mom1);
320       track2->GetPxPyPz(mom2);
321       //     Double_t cosTheta = (mom1[0]*mom2[0]+mom1[1]*mom2[1]+mom1[2]*mom2[2])/( TMath::Sqrt(mom1[0]*mom1[0]+mom1[1]*mom1[1]+mom1[2]*mom1[2])*TMath::Sqrt(mom2[0]*mom2[0]+mom2[1]*mom2[1]+mom2[2]*mom2[2]) );
322       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
323       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
324       //Double_t theta = momv1.Angle(momv2);
325       Double_t theta = momv1.Phi()-momv2.Phi();
326       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
327     
328       fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
329       Float_t deltaPhi = track1->Phi()-track2->Phi();
330       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
331       fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());
332
333       track1->GetImpactParameters(dcaR[0],dcaZ[0]);
334       track2->GetImpactParameters(dcaR[1],dcaZ[1]);
335
336       if(track2->Pt()<0.5) continue;
337       Double_t rIsol = TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) );
338       fRisol->Fill(rIsol); //Fill R histogram
339       if(track2->Pt()<fPtMin) continue;
340
341         fTheta->Fill(theta);
342         fThetaZoom->Fill(theta);
343         fThetaPt1Pt2->Fill(theta,track1->Pt(),track2->Pt());
344         fThetaPt1Pt2Signed->Fill(theta,track1->GetSign()*track1->Pt(),track2->GetSign()*track2->Pt());
345         fDeltaPhiSumEtaPt1->Fill(deltaPhi,track1->Eta()+track2->Eta(),track1->Pt());
346         fDeltaPhiSumEtaPt2->Fill(deltaPhi,track1->Eta()+track2->Eta(),track2->Pt());
347         fThetaDCAZ1DCAZ2->Fill(theta,dcaZ[0],dcaZ[1]);
348         fRisolTheta->Fill(rIsol,theta);
349         if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
350           fDCAZCosmicCandidates->Fill(dcaZ[0],dcaZ[1]);
351           fDCARCosmicCandidates->Fill(dcaR[0],dcaR[1]);
352           fPtSignedCosmicCandidates->Fill(track1->GetSign()*track1->Pt());
353         }
354
355     } // track2 loop
356     
357   } // track1 loop 
358
359   PostData(1,fHistListCosmics);
360 }      
361
362 //________________________________________________________________________
363 void AliPWG4CosmicCandidates::Terminate(Option_t *) 
364 {
365   //
366   // Called once at the end of the query
367   //
368   
369 }