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   Double_t *binsDCAR=new Double_t[fgkNDCARBins+1];
 192   for(Int_t i=0; i<=fgkNDCARBins; i++) binsDCAR[i]=(Double_t)fgkDCARMin + (fgkDCARMax-fgkDCARMin)/fgkNDCARBins*(Double_t)i ;
 193
 194   Int_t fgkNDCAZBins=80;
 195   Float_t fgkDCAZMin = -2.;
 196   Float_t fgkDCAZMax = 2.;
 197   Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
 198   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_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 }