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