]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliAnalysisTaskJetCore.cxx
Initial analysis task for TRD jet triggered data (J. Klein)
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
index 60ef95d672afa5c981770a855c694ec33749d307..d6ef8f9a8ac81df1cfd13d9c4ab9e2006bbb31bd 100644 (file)
@@ -74,8 +74,8 @@ fVtxZMax(10.),
 fEvtClassMin(0),
 fEvtClassMax(4),
 fFilterMask(0),
-fFilterMaskBestPt(16),
-fFilterType(2),
+fFilterMaskBestPt(0),
+fFilterType(0),
 fRadioFrac(0.2),
 fMinDist(0.1),
 fCentMin(0.),
@@ -90,10 +90,13 @@ fFlagEtaBkg(0),
 fFlagJetHadron(0),
 fFlagRandom(0),
 fFlagOnlyRecoil(0),
-fFlagOnlyHardest(0),
+fFlagOnlyHardest(1),
 fTrackTypeRec(kTrackUndef),
 fRPAngle(0),
-fNRPBins(100),
+fNRPBins(50),
+fSemigoodCorrect(0),
+fHolePos(4.71),
+fHoleWidth(0.2),
 fJetEtaMin(-.5),
 fJetEtaMax(.5),
 fNevents(0),
@@ -106,6 +109,7 @@ fJetPtFractionMin(0.5),
 fNMatchJets(4),
 fMatchMaxDist(0.8),
 fKeepJets(kFALSE),
+fRunAnaAzimuthalCorrelation(kFALSE),
 fkNbranches(2),
 fkEvtClasses(12),
 fOutputList(0x0),
@@ -140,14 +144,17 @@ fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),
 fh2Ntriggers(0x0),
-fh2Ntriggers2(0x0), 
+fh2Ntriggers2C10(0x0),
+fh2Ntriggers2C20(0x0), 
 fh3JetDensity(0x0),
 fh3JetDensityA4(0x0),
-fh2RPJets(0x0),
-fh2RPT(0x0),
-fh3spectriggeredC20RP(0x0),
-fh3spectriggeredC20(0x0),
-fh3spectriggeredC3060(0x0)
+fh2RPJetsC10(0x0),
+fh2RPJetsC20(0x0),
+fh2RPTC10(0x0),
+fh2RPTC20(0x0), 
+fHJetSpec(0x0),
+fhTTPt(0x0),
+fHJetPhiCorr(0x0)
 
  
 {
@@ -188,8 +195,8 @@ fVtxZMax(10.),
 fEvtClassMin(0),
 fEvtClassMax(4),
 fFilterMask(0),
-fFilterMaskBestPt(16),
-fFilterType(2),
+fFilterMaskBestPt(0),
+fFilterType(0),
 fRadioFrac(0.2),
 fMinDist(0.1),
 fCentMin(0.),
@@ -204,10 +211,13 @@ fFlagEtaBkg(0),
 fFlagJetHadron(0),
 fFlagRandom(0),
 fFlagOnlyRecoil(0),
-fFlagOnlyHardest(0),
+fFlagOnlyHardest(1),
 fTrackTypeRec(kTrackUndef),
 fRPAngle(0),
-fNRPBins(100),
+fNRPBins(50),
+fSemigoodCorrect(0),
+fHolePos(4.71),
+fHoleWidth(0.2),
 fJetEtaMin(-.5),
 fJetEtaMax(.5),
 fNevents(0),
@@ -220,6 +230,7 @@ fJetPtFractionMin(0.5),
 fNMatchJets(4),
 fMatchMaxDist(0.8),
 fKeepJets(kFALSE),
+fRunAnaAzimuthalCorrelation(kFALSE),
 fkNbranches(2),
 fkEvtClasses(12),
 fOutputList(0x0),
@@ -254,14 +265,17 @@ fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),    
 fh2Ntriggers(0x0),
-fh2Ntriggers2(0x0),
+fh2Ntriggers2C10(0x0),
+fh2Ntriggers2C20(0x0),
 fh3JetDensity(0x0),
 fh3JetDensityA4(0x0),
-fh2RPJets(0x0),
-fh2RPT(0x0),
-fh3spectriggeredC20RP(0x0),
-fh3spectriggeredC20(0x0),
-fh3spectriggeredC3060(0x0)
+fh2RPJetsC10(0x0),
+fh2RPJetsC20(0x0),
+fh2RPTC10(0x0),
+fh2RPTC20(0x0), 
+fHJetSpec(0x0),
+fhTTPt(0x0),
+fHJetPhiCorr(0x0)
 
  {
    // Constructor
@@ -396,15 +410,15 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
    
 
     fh2Ntriggers=new TH2F("# of triggers","",10,0.,100.,50,0.,50.);
-    fh2Ntriggers2=new TH2F("# of triggers2","",50,0.,50.,50,0.,50.);
-
+    fh2Ntriggers2C10=new TH2F("# of triggers2C10","",50,0.,50.,50,0.,50.);
+    fh2Ntriggers2C20=new TH2F("# of triggers2C20","",50,0.,50.,50,0.,50.);
     fh3JetDensity=new TH3F("Jet density vs mutliplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
     fh3JetDensityA4=new TH3F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
-    fh2RPJets=new TH2F("RPJet","",fNRPBins,0.,fNRPBins,150,0.,150.);
-    fh2RPT=new TH2F("RPTrigger","",fNRPBins,0.,fNRPBins,150,0.,150.); 
-    fh3spectriggeredC20RP = new TH3F("Triggered spectrumC20RP","",3,0.,3.,140,-80.,200.,10,0.,50.);
-    fh3spectriggeredC20 = new TH3F("Triggered spectrumC20","",100,0.,1.,140,-80.,200.,10,0.,50.);
-    fh3spectriggeredC3060 = new TH3F("Triggered spectrumC3060","",100,0.,1.,140,-80.,200.,10,0.,50.);
+    fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
+    fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.); 
+    fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.); 
+    fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,3.5,50,0.,50.);  
+
 
     
     
@@ -457,15 +471,38 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
 
  
        fOutputList->Add(fh2Ntriggers);
-        fOutputList->Add(fh2Ntriggers2);
+        fOutputList->Add(fh2Ntriggers2C10);
+        fOutputList->Add(fh2Ntriggers2C20); 
         fOutputList->Add(fh3JetDensity);
         fOutputList->Add(fh3JetDensityA4);
-        fOutputList->Add(fh2RPJets);
-        fOutputList->Add(fh2RPT);
-        fOutputList->Add(fh3spectriggeredC20RP); 
-       fOutputList->Add(fh3spectriggeredC20); 
-       fOutputList->Add(fh3spectriggeredC3060);   
+        fOutputList->Add(fh2RPJetsC10);
+        fOutputList->Add(fh2RPJetsC20);
+         fOutputList->Add(fh2RPTC10);
+        fOutputList->Add(fh2RPTC20);
+
+        const Int_t dimSpec = 5;
+       const Int_t nBinsSpec[dimSpec]     = {10,100, 140, 50, fNRPBins};
+       const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0, 0};
+       const Double_t hiBinSpec[dimSpec]  = {100,1, 200, 50, fNRPBins};
+       fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
+       fOutputList->Add(fHJetSpec);  
+
+
+       if(fRunAnaAzimuthalCorrelation)
+         {
+           fhTTPt = new TH2F("fhTTPt","Trigger track p_{T} vs centrality",10,0,100,100,0,100);
+           fOutputList->Add(fhTTPt);
+
+           const Int_t dimCor = 5;
+           const Int_t nBinsCor[dimCor]     = {50, 200, 100,              8,   10};
+           const Double_t lowBinCor[dimCor] = {0,  -50, -0.5*TMath::Pi(), 0,   0};
+           const Double_t hiBinCor[dimCor]  = {50, 150, 1.5*TMath::Pi(),  0.8, 100};
+           fHJetPhiCorr = new THnSparseF("fHJetPhiCorr","TT p_{T} vs jet p_{T} vs dPhi vs area vs centrality",dimCor,nBinsCor,lowBinCor,hiBinCor);
+           fOutputList->Add(fHJetPhiCorr);
+         }
 
+
+        
      
    // =========== Switch on Sumw2 for all histos ===========
    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
@@ -658,30 +695,40 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    Int_t trigInTrack=-1;
    fRPAngle = aod->GetHeader()->GetEventplane();     
 
-   //   AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);     
-   //if(!partback){  
-   //PostData(1, fOutputList);
-   //return;}
-
-
-   for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
+   AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);     
+   if(!partback){  
+   PostData(1, fOutputList);
+   return;}
 
-   if(fFlagOnlyHardest!=0){if(tt!=nT) continue;}
-   AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);     
-   fh2Ntriggers->Fill(centValue,partback->Pt());
-   Int_t phiBinT = GetPhiBin(RelativePhi(partback->Phi(),fRPAngle));
-   if(centValue<20.) fh2RPT->Fill(phiBinT,partback->Pt()); 
 
+   //for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
+   //if(fFlagOnlyHardest!=0){if(tt!=nT) continue;}
+   //AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);     
    Double_t accep=2.*TMath::Pi()*1.8;
    Int_t injet4=0;
    Int_t injet=0; 
+   if(fSemigoodCorrect){
+   Double_t disthole=RelativePhi(partback->Phi(),fHolePos);
+   if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-0.6){ 
+   PostData(1, fOutputList);
+   return;}
+
+   }
+
+   fh2Ntriggers->Fill(centValue,partback->Pt());
+   Double_t phiBinT = RelativePhi(partback->Phi(),fRPAngle);
+   if(centValue<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
+   if(centValue<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
+
+
+
    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
            etabig  = jetbig->Eta();
            phibig  = jetbig->Phi();
            ptbig   = jetbig->Pt();
            if(ptbig==0) continue; 
-           Int_t phiBin = GetPhiBin(RelativePhi(phibig,fRPAngle));       
+           Double_t phiBin = RelativePhi(phibig,fRPAngle);       
            areabig = jetbig->EffectiveAreaCharged();
            Double_t ptcorr=ptbig-rho*areabig;
           if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
@@ -707,33 +754,40 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
            if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
 
 
-                   if(centValue<20.) fh2RPJets->Fill(phiBin, ptcorr);
-                
+          if(centValue<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
+          if(centValue<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
                    Double_t dismin=100.;
                    Double_t ptmax=-10.; 
                    Int_t index1=-1;
                    Int_t index2=-1;
+         
+                  Float_t phitt=partback->Phi();
+                   if(phitt<0)phitt+=TMath::Pi()*2.; 
+                   Int_t phiBintt = GetPhiBin(phitt-fRPAngle);
+
+                  Double_t fillspec[] = {centValue,jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt(),phiBintt};
+                 fHJetSpec->Fill(fillspec);
+           
           
-           if(centValue<20.)  fh3spectriggeredC20RP->Fill(phiBin,ptcorr,partback->Pt());
-           if(centValue<20.)  fh3spectriggeredC20->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
-           if(centValue>30. && centValue<60.)  fh3spectriggeredC3060->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
 
                    if(ptcorr<=0) continue;
 
-                  
                        AliAODTrack* leadtrack=0; 
                        Int_t ippt=0;
                        Double_t ppt=-10;
                        if(fFlagJetHadron==0){   
-                      TRefArray *genTrackList = jetbig->GetRefTracks();
-                       Int_t nTracksGenJet = genTrackList->GetEntriesFast();
-                       AliAODTrack* genTrack;
-                       for(Int_t ir=0; ir<nTracksGenJet; ++ir){
-                       genTrack = (AliAODTrack*)(genTrackList->At(ir));
-                      if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
-                      ippt=ir;}}
-                       leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
-                       if(!leadtrack) continue;}
+                        TRefArray *genTrackList = jetbig->GetRefTracks();
+                        Int_t nTracksGenJet = genTrackList->GetEntriesFast();
+                        AliAODTrack* genTrack;
+                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
+                          genTrack = (AliAODTrack*)(genTrackList->At(ir));
+                          if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
+                            ippt=ir;}}
+                        leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
+                        if(!leadtrack) continue;
+                      }
+
+
 
                       AliVParticle* leadtrackb=0;
                        if(fFlagJetHadron!=0){
@@ -746,26 +800,26 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
 
 
                        
-                       //store one trigger info                   
-                        if(iCount==0){                        
-                         trigJet=i;
-                          trigBBTrack=nT;
-                          trigInTrack=ippt;
-                          iCount=iCount+1;} 
-
+                      //store one trigger info                   
+                      if(iCount==0){                        
+                        trigJet=i;
+                        trigBBTrack=nT;
+                        trigInTrack=ippt;
+                        iCount=iCount+1;} 
+                      
    
-                 if(fCheckMethods){
-                  for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
-                  AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
-                  etasmall  = jetsmall->Eta();
-                  phismall = jetsmall->Phi();
-                  ptsmall   = jetsmall->Pt();
-                  areasmall = jetsmall->EffectiveAreaCharged();
-                  Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
-                 tmpDeltaR=TMath::Sqrt(tmpDeltaR);
-                    //Fraction in the jet core  
-                    if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
-                   index2=j;}  
+                      if(fCheckMethods){
+                        for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
+                          AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
+                          etasmall  = jetsmall->Eta();
+                          phismall = jetsmall->Phi();
+                          ptsmall   = jetsmall->Pt();
+                          areasmall = jetsmall->EffectiveAreaCharged();
+                          Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
+                          tmpDeltaR=TMath::Sqrt(tmpDeltaR);
+                          //Fraction in the jet core  
+                          if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
+                            index2=j;}  
                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
                      index1=j;}} //en of loop over R=0.2 jets
                 //method1:most concentric jet=core 
@@ -781,26 +835,25 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig); 
                  if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
                  if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}  
-       
-                  
-        
-        // for(int it = 0;it<ParticleList.GetEntries();++it){
-        //  AliVParticle *part = (AliVParticle*)ParticleList.At(it);
-                //  Double_t deltaR = jetbig->DeltaR(part);
-         //  Double_t deltaEta = etabig-part->Eta();
-         //  if(centValue<20.) fh2Ntriggers2->Fill(partback->Pt(),part->Pt());
-         //  if(fDoEventMixing==0 && fFlagOnlyRecoil==0){ 
-         //  Double_t deltaPhi=phibig-part->Phi();
-         //  if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
-         //  if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
-        //  Double_t pTcont=0;
-         //  if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
-         //  if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt(); 
-        //   Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};  
-         //   fhnDeltaR->Fill(jetEntries);}
-
-
-        // }
+                 if(centValue<10&&leadtrack) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());                  
+                  if(centValue<20&&leadtrack) fh2Ntriggers2C20->Fill(leadtrack->Pt(),partback->Pt());  
+         if(fDoEventMixing==0 && fFlagOnlyRecoil==0){ 
+        for(int it = 0;it<ParticleList.GetEntries();++it){
+         AliVParticle *part = (AliVParticle*)ParticleList.At(it);
+                 Double_t deltaR = jetbig->DeltaR(part);
+          Double_t deltaEta = etabig-part->Eta();
+         
+          Double_t deltaPhi=phibig-part->Phi();
+          if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
+          if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
+         Double_t pTcont=0;
+          if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
+          if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt(); 
+          Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};  
+           fhnDeltaR->Fill(jetEntries);}
+
+
+         }
         
         
           //end of track loop, we only do it if EM is switched off
@@ -818,7 +871,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    if(injet4>0)fh3JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep,partback->Pt());
           //end of jet loop
 
-   }
+   //}
 
 
           if(fDoEventMixing>0){
@@ -869,6 +922,33 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
          
          }
 
+         /////////////////////////////////////////////////////////////////////////////
+         ////////////////////// Rongrong's analysis //////////////////////////////////
+         if(fRunAnaAzimuthalCorrelation)
+           {
+             fhTTPt->Fill(centValue,partback->Pt());
+             for(Int_t ij=0; ij<fListJets[0]->GetEntries(); ij++)
+               {
+                 AliAODJet* jet = (AliAODJet*)(fListJets[0]->At(ij));
+                 Double_t jetPt   = jet->Pt();
+                 Double_t jetEta  = jet->Eta();
+                 Double_t jetPhi  = jet->Phi();
+                 if(jetPt==0) continue; 
+                 if((jetEta<fJetEtaMin)||(jetEta>fJetEtaMax)) continue;
+                 Double_t jetArea = jet->EffectiveAreaCharged();
+                 Double_t jetPtCorr=jetPt-rho*jetArea;
+                 Double_t dPhi=jetPhi-partback->Phi();
+                 if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
+                 if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
+                 if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
+                 if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
+
+                 Double_t fill[] = {partback->Pt(),jetPtCorr,dPhi,jetArea,centValue};
+                 fHJetPhiCorr->Fill(fill);
+               }
+           }
+         /////////////////////////////////////////////////////////////////////////////
+         /////////////////////////////////////////////////////////////////////////////
 
 
      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
@@ -990,7 +1070,10 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
 
      if(!fESD)aod = fAODIn;
      else aod = fAODOut;   
-      Int_t index=-1;
+
+     if(!aod)return 0;
+
+     Int_t index=-1;
      Double_t ptmax=-10;
 
 
@@ -1000,14 +1083,14 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
       Bool_t bGood = false;
       if(fFilterType == 0)bGood = true;
       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
-      else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();       
-      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
+      else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();    
+     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
       if(bGood==false) continue;
       if(TMath::Abs(tr->Eta())>0.9)continue;
       if(tr->Pt()<0.15)continue;
       list->Add(tr);
       iCount++;
-      if(fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
+      if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
        if(tr->TestFilterBit(fFilterMaskBestPt)){
          if(tr->Pt()>ptmax){ 
            ptmax=tr->Pt();