Charged jets (pPb): Convert static to class member
[u/mrichter/AliRoot.git] / PWGJE / 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   Int_t fgkNPtBins2D= (int)((fgkPtMax-fgkPtMin)/4.);
185
186   Int_t fgkNPhiBins=18;
187   Float_t kMinPhi = -0.5*TMath::Pi();
188   Float_t kMaxPhi = 3./2.*TMath::Pi();
189
190   Int_t fgkNThetaBins=fgkNPhiBins*8;
191   Float_t kMinTheta = -0.5*TMath::Pi();
192   Float_t kMaxTheta = 3./2.*TMath::Pi();
193
194   Int_t fgkNDCARBins=40;
195   Float_t fgkDCARMin = -0.2;
196   Float_t fgkDCARMax = 0.2;
197   Float_t *binsDCAR=new Float_t[fgkNDCARBins+1];
198   for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Float_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Float_t)i ;
199
200   Int_t fgkNDCAZBins=40;
201   Float_t fgkDCAZMin = -2.;
202   Float_t fgkDCAZMax = 2.;
203   Float_t *binsDCAZ=new Float_t[fgkNDCAZBins+1];
204   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Float_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Float_t)i ;
205
206   fPtSignedCosmicCandidates = new TH1F("fPtSignedCosmicCandidates","fPtSignedCosmicCandidates",2*fgkNPtBins, -1.*fgkPtMax, fgkPtMax);
207   fHistListCosmics->Add(fPtSignedCosmicCandidates);  
208
209   fDeltaPtCosmicCandidates = new TH1F("fDeltaPtCosmicCandidates","fDeltaPtCosmicCandidates",fgkNPtBins, -50., 50.);
210   fHistListCosmics->Add(fDeltaPtCosmicCandidates);  
211
212   fDeltaPhiSumEta = new TH2F("fDeltaPhiSumEta","fDeltaPhiSumEta",fgkNPhiBins*4,kMinPhi,kMaxPhi,80, -2.,2.);
213   fHistListCosmics->Add(fDeltaPhiSumEta);  
214
215   fDCAZCosmicCandidates = new TH2F("fDCAZCosmicCandidates","fDCAZCosmicCandidates",fgkNDCAZBins,binsDCAZ,fgkNDCAZBins,binsDCAZ);
216   fHistListCosmics->Add(fDCAZCosmicCandidates);
217
218   fDCARCosmicCandidates = new TH2F("fDCARCosmicCandidates","fDCARCosmicCandidates",fgkNDCARBins,binsDCAR,fgkNDCARBins,binsDCAR);
219   fHistListCosmics->Add(fDCARCosmicCandidates);
220
221   fTheta = new TH1F("fTheta","fTheta",fgkNThetaBins,kMinTheta,kMaxTheta);
222   fHistListCosmics->Add(fTheta);
223
224   fThetaZoom = new TH1F("fThetaZoom","fThetaZoom",100,TMath::Pi()-1.,TMath::Pi()+1.);
225   fHistListCosmics->Add(fThetaZoom);
226
227   fThetaPt1Pt2 = new TH3F("fThetaPt1Pt2","fThetaPt1Pt2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,fgkPtMin,fgkPtMax,fgkNPtBins2D,fgkPtMin,fgkPtMax);
228   fHistListCosmics->Add(fThetaPt1Pt2);
229
230   fThetaPt1Pt2Signed = new TH3F("fThetaPt1Pt2Signed","fThetaPt1Pt2Signed",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNPtBins2D,-1.*fgkPtMax,fgkPtMax,fgkNPtBins2D,-1.*fgkPtMax,fgkPtMax);
231   fHistListCosmics->Add(fThetaPt1Pt2Signed);
232
233   fDeltaPhiSumEtaPt1 = new TH3F("fDeltaPhiSumEtaPt1","fDeltaPhiSumEtaPt1",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
234   fHistListCosmics->Add(fDeltaPhiSumEtaPt1);
235
236   fDeltaPhiSumEtaPt2 = new TH3F("fDeltaPhiSumEtaPt2","fDeltaPhiSumEtaPt2",fgkNThetaBins,kMinTheta,kMaxTheta,80, -2.,2.,fgkNPtBins2D,fgkPtMin,fgkPtMax);
237   fHistListCosmics->Add(fDeltaPhiSumEtaPt2);
238
239   fThetaDCAZ1DCAZ2 = new TH3F("fThetaDCAZ1DCAZ2","fThetaDCAZ1DCAZ2",fgkNThetaBins,kMinTheta,kMaxTheta,fgkNDCAZBins,-2.,2.,fgkNDCAZBins,-2.,2.);
240   fHistListCosmics->Add(fThetaDCAZ1DCAZ2);
241
242   fRisol = new TH1F("fRisol","fRisol",100,0.,10.);
243   fHistListCosmics->Add(fRisol);
244
245   fRisolTheta = new TH2F("fRisolTheta","fRisolTheta",100,0.,10.,fgkNThetaBins,kMinTheta,kMaxTheta);
246   fHistListCosmics->Add(fRisolTheta);
247
248   TH1::AddDirectory(oldStatus); 
249
250   PostData(1, fHistListCosmics);  
251
252   if(binsDCAR) delete [] binsDCAR;
253   if(binsDCAZ) delete [] binsDCAZ;
254
255 }
256
257 //________________________________________________________________________
258 void AliPWG4CosmicCandidates::UserExec(Option_t *) 
259 {
260   // Main loop
261   // Called for each event
262   //
263
264   // All events without selection
265   fNEventAll->Fill(0.);
266    
267   if (!fInputEvent) {
268     AliDebug(2,Form("ERROR: fESD not available"));
269     cout << "ERROR: fESD not available" << endl;
270     PostData(1, fHistListCosmics);
271     return;
272   }
273
274   const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
275   if(!vtx){
276     PostData(1, fHistListCosmics);
277     return;
278   }
279
280 // Need vertex cut
281   TString vtxName(vtx->GetName());
282   if( vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
283     // Post output data
284     PostData(1, fHistListCosmics);
285     return;
286   }
287
288   //  AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
289   double primVtx[3];
290   vtx->GetXYZ(primVtx);
291   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
292     // Post output data
293     PostData(1, fHistListCosmics);
294     return;
295   }
296   if(!fInputEvent->GetNumberOfTracks() || fInputEvent->GetNumberOfTracks()<2){ 
297     // Post output data
298     PostData(1, fHistListCosmics);
299     return;
300   }
301   Int_t nTracks = fInputEvent->GetNumberOfTracks();
302
303   if(!fTrackCuts) {
304    // Post output data
305     PostData(1, fHistListCosmics);
306     return;
307   }
308
309   fNEventSel->Fill(0.);
310
311   Float_t dcaR[2] = {0.,0.};
312   Float_t dcaZ[2] = {0.,0.};
313
314   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
315
316     AliESDtrack* track1 = (AliESDtrack*)fInputEvent->GetTrack(iTrack1);
317     if (!track1)  continue;
318     if(!(fTrackCuts->AcceptTrack(track1))) { continue; }
319     if(track1->Pt()<fPtMin) continue;
320     //Start 2nd track loop to look for correlations
321     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
322       AliESDtrack *track2 = (AliESDtrack*)fInputEvent->GetTrack(iTrack2);
323       if(!track2) continue;
324       if(!(fTrackCuts->AcceptTrack(track2))) { continue; }
325          
326       //Check if back-to-back
327       Double_t mom1[3],mom2[3];
328       track1->GetPxPyPz(mom1);
329       track2->GetPxPyPz(mom2);
330       //     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]) );
331       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
332       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
333       //Double_t theta = momv1.Angle(momv2);
334       Double_t theta = momv1.Phi()-momv2.Phi();
335       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
336     
337       fDeltaPtCosmicCandidates->Fill(track1->Pt()-track2->Pt());
338       Float_t deltaPhi = track1->Phi()-track2->Phi();
339       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
340       fDeltaPhiSumEta->Fill(deltaPhi,track1->Eta()+track2->Eta());
341
342       track1->GetImpactParameters(dcaR[0],dcaZ[0]);
343       track2->GetImpactParameters(dcaR[1],dcaZ[1]);
344
345       if(track2->Pt()<0.5) continue;
346       Double_t rIsol = TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) );
347       fRisol->Fill(rIsol); //Fill R histogram
348       if(track2->Pt()<fPtMin) continue;
349
350         fTheta->Fill(theta);
351         fThetaZoom->Fill(theta);
352         fThetaPt1Pt2->Fill(theta,track1->Pt(),track2->Pt());
353         fThetaPt1Pt2Signed->Fill(theta,track1->GetSign()*track1->Pt(),track2->GetSign()*track2->Pt());
354         fDeltaPhiSumEtaPt1->Fill(deltaPhi,track1->Eta()+track2->Eta(),track1->Pt());
355         fDeltaPhiSumEtaPt2->Fill(deltaPhi,track1->Eta()+track2->Eta(),track2->Pt());
356         fThetaDCAZ1DCAZ2->Fill(theta,dcaZ[0],dcaZ[1]);
357         fRisolTheta->Fill(rIsol,theta);
358         if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
359           fDCAZCosmicCandidates->Fill(dcaZ[0],dcaZ[1]);
360           fDCARCosmicCandidates->Fill(dcaR[0],dcaR[1]);
361           fPtSignedCosmicCandidates->Fill(track1->GetSign()*track1->Pt());
362         }
363
364     } // track2 loop
365     
366   } // track1 loop 
367
368   PostData(1,fHistListCosmics);
369 }      
370
371 //________________________________________________________________________
372 void AliPWG4CosmicCandidates::Terminate(Option_t *) 
373 {
374   //
375   // Called once at the end of the query
376   //
377   
378 }