b17fad059da4722482a7d90bb7efb537e6c7b811
[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 //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     // Dummy copy constructor
140 }
141
142 //________________________________________________________________________
143 AliPWG4CosmicCandidates& AliPWG4CosmicCandidates::operator=(const AliPWG4CosmicCandidates& /*trclass*/)
144 {
145     // Dummy assignment operator
146     return *this;
147 }
148
149 //________________________________________________________________________
150 void AliPWG4CosmicCandidates::LocalInit()
151 {
152   //
153   // Only called once at beginning
154   //
155   PostData(2,fTrackCuts);
156 }
157
158 //________________________________________________________________________
159 void AliPWG4CosmicCandidates::UserCreateOutputObjects()
160 {
161   //Create output objects
162   // Called once
163   AliDebug(2,Form(">> AliPWG4CosmicCandidates::UserCreateOutputObjects \n")); 
164
165   Bool_t oldStatus = TH1::AddDirectoryStatus();
166   TH1::AddDirectory(kFALSE); 
167   
168   OpenFile(1);
169   fHistListCosmics = new TList();
170
171   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
172   fHistListCosmics->Add(fNEventAll);
173   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
174   fHistListCosmics->Add(fNEventSel);
175
176   Float_t fgkPtMin=0.;
177   Float_t fgkPtMax=100.;
178   Int_t fgkNPtBins= (int)(fgkPtMax-fgkPtMin);
179
180   Int_t fgkNPhiBins=18;
181   Float_t kMinPhi = -0.5*TMath::Pi();
182   Float_t kMaxPhi = 3./2.*TMath::Pi();
183
184   Int_t fgkNThetaBins=fgkNPhiBins*8;
185   Float_t kMinTheta = -0.5*TMath::Pi();
186   Float_t kMaxTheta = 3./2.*TMath::Pi();
187
188   Int_t fgkNDCARBins=80;
189   Float_t fgkDCARMin = -0.2;
190   Float_t fgkDCARMax = 0.2;
191   Float_t *binsDCAR=new Float_t[fgkNDCARBins+1];
192   for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Float_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Float_t)i ;
193
194   Int_t fgkNDCAZBins=80;
195   Float_t fgkDCAZMin = -2.;
196   Float_t fgkDCAZMax = 2.;
197   Float_t *binsDCAZ=new Float_t[fgkNDCAZBins+1];
198   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Float_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Float_t)i ;
199
200   fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2*(int)(fgkPtMax-fgkPtMin), -1.*fgkPtMax, fgkPtMax);
201   fHistListCosmics->Add(fPtSignedCosmicCandidates);  
202
203   fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
204   fHistListCosmics->Add(fDeltaPtCosmicCandidates);  
205
206   fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,kMinPhi,kMaxPhi,80, -2.,2.);
207   fHistListCosmics->Add(fDeltaPhiSumEta);  
208
209   fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
210   fHistListCosmics->Add(fDCAZCosmicCandidates);
211
212   fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
213   fHistListCosmics->Add(fDCARCosmicCandidates);
214
215   fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
216   fHistListCosmics->Add(fTheta);
217
218   fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
219   fHistListCosmics->Add(fThetaZoom);
220
221   fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
222   fHistListCosmics->Add(fThetaPt1Pt2);
223
224   fThetaPt1Pt2Signed = new TH3F("fThetaPt1Pt2Signed","fThetaPt1Pt2Signed",fgkNThetaBins,kMinTheta,kMaxTheta,4*(int)(fgkPtMax-fgkPtMin),-1.*fgkPtMax,fgkPtMax,4*(int)(fgkPtMax-fgkPtMin),-1.*fgkPtMax,fgkPtMax);
225   fHistListCosmics->Add(fThetaPt1Pt2Signed);
226
227   fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
228   fHistListCosmics->Add(fDeltaPhiSumEtaPt1);
229
230   fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,(int)(fgkPtMax-fgkPtMin),fgkPtMin,fgkPtMax);
231   fHistListCosmics->Add(fDeltaPhiSumEtaPt2);
232
233   fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
234   fHistListCosmics->Add(fThetaDCAZ1DCAZ2);
235
236   fRisol = new TH1F("fRisol","fRisol",100,0.,10.);
237   fHistListCosmics->Add(fRisol);
238
239   fRisolTheta = new TH2F("fRisolTheta","fRisolTheta",100,0.,10.,fgkNThetaBins,kMinTheta,kMaxTheta);
240   fHistListCosmics->Add(fRisolTheta);
241
242   TH1::AddDirectory(oldStatus);   
243
244   if(binsDCAR) delete [] binsDCAR;
245   if(binsDCAZ) delete [] binsDCAZ;
246
247 }
248
249 //________________________________________________________________________
250 void AliPWG4CosmicCandidates::UserExec(Option_t *) 
251 {
252 //   // Main loop
253 //   // Called for each event
254
255   // All events without selection
256   fNEventAll->Fill(0.);
257    
258   if (!fInputEvent) {
259     AliDebug(2,Form("ERROR: fESD not available"));
260     cout << "ERROR: fESD not available" << endl;
261     PostData(1, fHistListCosmics);
262     return;
263   }
264
265   const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
266   // Need vertex cut
267   if (!vtx || vtx->GetNContributors() < 2) {
268     // Post output data
269     PostData(1, fHistListCosmics);
270     return;
271   }
272
273   //  AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
274   double primVtx[3];
275   vtx->GetXYZ(primVtx);
276   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
277     // Post output data
278     PostData(1, fHistListCosmics);
279     return;
280   }
281   if(!fInputEvent->GetNumberOfTracks() || fInputEvent->GetNumberOfTracks()<2){ 
282     // Post output data
283     PostData(1, fHistListCosmics);
284     return;
285   }
286   Int_t nTracks = fInputEvent->GetNumberOfTracks();
287
288   if(!fTrackCuts) {
289    // Post output data
290     PostData(1, fHistListCosmics);
291     return;
292   }
293
294  fNEventSel->Fill(0.);
295
296   Float_t dcaR[2] = {0.,0.};
297   Float_t dcaZ[2] = {0.,0.};
298
299   // Track loop to fill a pT spectrum
300   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
301
302     AliESDtrack* track1 = (AliESDtrack*)fInputEvent->GetTrack(iTrack1);
303     if (!track1)  continue;
304     if(!(fTrackCuts->AcceptTrack(track1))) { continue; }
305     if(track1->Pt()<fPtMin) continue;
306     //Start 2nd track loop to look for correlations
307     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
308       AliESDtrack *track2 = (AliESDtrack*)fInputEvent->GetTrack(iTrack2);
309       if(!track2) continue;
310       if(!(fTrackCuts->AcceptTrack(track2))) { continue; }
311          
312       //Check if back-to-back
313       Double_t mom1[3],mom2[3];
314       track1->GetPxPyPz(mom1);
315       track2->GetPxPyPz(mom2);
316       //     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]) );
317       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
318       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
319       //Double_t theta = momv1.Angle(momv2);
320       Double_t theta = momv1.Phi()-momv2.Phi();
321       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
322     
323       fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
324       Float_t deltaPhi = track1->Phi()-track2->Phi();
325       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
326       fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());
327
328       track1->GetImpactParameters(dcaR[0],dcaZ[0]);
329       track2->GetImpactParameters(dcaR[1],dcaZ[1]);
330
331       if(track2->Pt()<0.5) continue;
332       Double_t rIsol = TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) );
333       fRisol->Fill(rIsol); //Fill R histogram
334       if(track2->Pt()<fPtMin) continue;
335
336         fTheta->Fill(theta);
337         fThetaZoom->Fill(theta);
338         fThetaPt1Pt2->Fill(theta,track1->Pt(),track2->Pt());
339         fThetaPt1Pt2Signed->Fill(theta,track1->GetSign()*track1->Pt(),track2->GetSign()*track2->Pt());
340         fDeltaPhiSumEtaPt1->Fill(deltaPhi,track1->Eta()+track2->Eta(),track1->Pt());
341         fDeltaPhiSumEtaPt2->Fill(deltaPhi,track1->Eta()+track2->Eta(),track2->Pt());
342         fThetaDCAZ1DCAZ2->Fill(theta,dcaZ[0],dcaZ[1]);
343         fRisolTheta->Fill(rIsol,theta);
344         if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
345           fDCAZCosmicCandidates->Fill(dcaZ[0],dcaZ[1]);
346           fDCARCosmicCandidates->Fill(dcaR[0],dcaR[1]);
347           fPtSignedCosmicCandidates->Fill(track1->GetSign()*track1->Pt());
348         }
349
350     } // track2 loop
351     
352   } // track1 loop 
353
354   PostData(1,fHistListCosmics);
355 }      
356
357 //________________________________________________________________________
358 void AliPWG4CosmicCandidates::Terminate(Option_t *) 
359 {
360   // Called once at the end of the query
361
362   
363 }