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 }