Add GeneratorOnly option and hard-code asymmetric dE/dx cut
authordgangadh <dhevan.raja.gangadharan@cern.ch>
Mon, 13 Jan 2014 12:45:45 +0000 (13:45 +0100)
committerdgangadh <dhevan.raja.gangadharan@cern.ch>
Mon, 13 Jan 2014 12:45:45 +0000 (13:45 +0100)
PWGCF/FEMTOSCOPY/Chaoticity/AliChaoticity.cxx
PWGCF/FEMTOSCOPY/Chaoticity/AliChaoticity.h

index f6f5510..ecbe8d1 100644 (file)
@@ -57,6 +57,7 @@ AliAnalysisTaskSE(),
   fAODcase(kTRUE),
   fPbPbcase(kTRUE),
   fGenerateSignal(kFALSE),
+  fGeneratorOnly(kFALSE),
   fPdensityExplicitLoop(kFALSE),
   fPdensityPairCut(kTRUE),
   fTabulatePairs(kFALSE),
@@ -220,6 +221,7 @@ AliChaoticity::AliChaoticity(const Char_t *name)
   fAODcase(kTRUE),
   fPbPbcase(kTRUE),
   fGenerateSignal(kFALSE),
+  fGeneratorOnly(kFALSE),
   fPdensityExplicitLoop(kFALSE),
   fPdensityPairCut(kTRUE),
   fTabulatePairs(kFALSE),
@@ -389,6 +391,7 @@ AliChaoticity::AliChaoticity(const AliChaoticity &obj)
     fAODcase(obj.fAODcase),
     fPbPbcase(obj.fPbPbcase),
     fGenerateSignal(obj.fGenerateSignal),
+    fGeneratorOnly(obj.fGeneratorOnly),
     fPdensityExplicitLoop(obj.fPdensityExplicitLoop),
     fPdensityPairCut(obj.fPdensityPairCut),
     fTabulatePairs(obj.fTabulatePairs),
@@ -503,6 +506,7 @@ AliChaoticity &AliChaoticity::operator=(const AliChaoticity &obj)
   fAODcase = obj.fAODcase;
   fPbPbcase = obj.fPbPbcase; 
   fGenerateSignal = obj.fGenerateSignal;
+  fGeneratorOnly = obj.fGeneratorOnly;
   fPdensityExplicitLoop = obj.fPdensityExplicitLoop;
   fPdensityPairCut = obj.fPdensityPairCut;
   fTabulatePairs = obj.fTabulatePairs;
@@ -777,7 +781,7 @@ void AliChaoticity::ParInit()
   fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeights);
   for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
   //
-  fDampStart = 0.5;// was 0.3
+  fDampStart = 0.5;// was 0.3, then 0.5
   fDampStep = 0.02;
   
   //
@@ -1581,17 +1585,17 @@ void AliChaoticity::Exec(Option_t *)
       Double_t integratedTimesTOF[10]={0};
 
       /*if(fFilterBit != 7 ) {
-       nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron));
-       nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon));
-       nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion));
-       nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon));
-       nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton));
+       nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron);
+       nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon);
+       nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion);
+       nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon);
+       nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton);
        //
-       nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron));
-       nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon));
-       nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion));
-       nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon));
-       nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton));
+       nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron);
+       nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon);
+       nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion);
+       nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon);
+       nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton);
        signalTPC = aodtrack->GetTPCsignal();
        if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
          fTempStruct[myTracks].fTOFhit = kTRUE;
@@ -1608,17 +1612,17 @@ void AliChaoticity::Exec(Option_t *)
          
          UInt_t status2=aodTrack2->GetStatus();
          
-         nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron));
-         nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon));
-         nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion));
-         nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon));
-         nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton));
+         nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron);
+         nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon);
+         nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion);
+         nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon);
+         nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton);
          //
-         nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron));
-         nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon));
-         nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion));
-         nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon));
-         nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton));
+         nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron);
+         nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon);
+         nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion);
+         nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon);
+         nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton);
          signalTPC = aodTrack2->GetTPCsignal();
          
          if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
@@ -1640,15 +1644,17 @@ void AliChaoticity::Exec(Option_t *)
       
       // Use TOF if good hit and above threshold
       if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
-       if(nSigmaTOF[0]<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
-       if(nSigmaTOF[2]<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
-       if(nSigmaTOF[3]<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
-       if(nSigmaTOF[4]<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
+       if(fabs(nSigmaTOF[0])<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
+       if(fabs(nSigmaTOF[2])<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
+       if(fabs(nSigmaTOF[3])<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
+       if(fabs(nSigmaTOF[4])<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
       }else {// TPC info instead
-       if(nSigmaTPC[0]<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
-       if(nSigmaTPC[2]<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
-       if(nSigmaTPC[3]<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
-       if(nSigmaTPC[4]<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
+       if(fabs(nSigmaTPC[0])<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
+       if(fabs(nSigmaTPC[2])<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
+       if(fabs(nSigmaTPC[3])<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
+       if(fabs(nSigmaTPC[4])<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
+       // asymmetric dE/dx cut to reduce muon contamination
+       if(nSigmaTPC[2] < -0.5) fTempStruct[myTracks].fPion = kFALSE;// asymmetric dE/dx cut
       }
                
       
@@ -1711,6 +1717,46 @@ void AliChaoticity::Exec(Option_t *)
     return;
   }
   
+  // Generator info only
+  if(fMCcase && fGeneratorOnly){
+    myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
+    for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
+      if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
+      if(myTracks >= 1300) continue;// additional cut to limit high mult events which exceed pair # limits
+      
+      AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
+      if(!mcParticle) continue;
+      if(fabs(mcParticle->Eta())>0.8) continue;
+      if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
+      if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
+      if(!mcParticle->IsPrimary()) continue;
+      if(!mcParticle->IsPhysicalPrimary()) continue;
+      if(abs(mcParticle->GetPdgCode())!=211) continue;
+      
+      fTempStruct[myTracks].fP[0] = mcParticle->Px();
+      fTempStruct[myTracks].fP[1] = mcParticle->Py();
+      fTempStruct[myTracks].fP[2] = mcParticle->Pz();
+      fTempStruct[myTracks].fX[0] = 0.; fTempStruct[myTracks].fX[1] = 0.; fTempStruct[myTracks].fX[2] = 0.;
+      
+      fTempStruct[myTracks].fId = myTracks;// use my track counter 
+      fTempStruct[myTracks].fLabel = mctrackN;
+      fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
+      if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
+      fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
+      fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
+      fTempStruct[myTracks].fEta = mcParticle->Eta();
+      fTempStruct[myTracks].fCharge = int(mcParticle->Charge()/3.);
+      fTempStruct[myTracks].fDCAXY = 0.;
+      fTempStruct[myTracks].fDCAZ = 0.;
+      fTempStruct[myTracks].fDCA = 0.;
+      fTempStruct[myTracks].fPion = kTRUE;
+      fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2)); 
+      fTempStruct[myTracks].fKey = 1;
+      
+      myTracks++;
+      pionCount++;
+    }
+  }
   
   if(myTracks >= 1) {
     ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
@@ -1973,7 +2019,7 @@ void AliChaoticity::Exec(Option_t *)
        
        // Pair Splitting/Merging cut
        if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
-       if(ch1 == ch2){
+       if(ch1 == ch2 && !fGeneratorOnly){
          if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
            fPairSplitCut[0][i]->AddAt('1',j);
            ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
@@ -1998,7 +2044,7 @@ void AliChaoticity::Exec(Option_t *)
              ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
              ((TH1D*)fOutputList->FindObject("fQDist"))->Fill(qinv12-qinv12MC);
            }
-           
+           //cout<<pVect1[0]<<"  "<<pVect1MC[0]<<" : "<<pVect1[1]<<"  "<<pVect1MC[1]<<" : "<<pVect1[2]<<"  "<<pVect1MC[2]<<" : "<<pVect1[3]<<"  "<<pVect1MC[3]<<" : "<<endl;
            //if(transK12 <= 0.35) fEDbin=0;
            //else fEDbin=1;
 
@@ -2085,6 +2131,7 @@ void AliChaoticity::Exec(Option_t *)
            if(fGenerateSignal) {
              WInput = MCWeight(ch1,ch2, fRMax, fFixedLambdaBinr3, qinv12, transK12);
              //WInput = MCWeight(ch1,ch2, fRMax, fFixedLambdaBinr3, qinv12MC);
+             //cout<<qinv12<<"  "<<qinv12MC<<endl;
              KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong), WInput);
            }else KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
            
@@ -2492,7 +2539,7 @@ void AliChaoticity::Exec(Option_t *)
        
 
        if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
-       if(ch1 == ch2){
+       if(ch1 == ch2 && !fGeneratorOnly){
          if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
            fPairSplitCut[1][i]->AddAt('1',j);
            continue;
@@ -2658,7 +2705,7 @@ void AliChaoticity::Exec(Option_t *)
        ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
        ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
        
-       if(ch1 == ch2){
+       if(ch1 == ch2 && !fGeneratorOnly){
          if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
            fPairSplitCut[2][i]->AddAt('1',j);
            continue;
@@ -2711,7 +2758,7 @@ void AliChaoticity::Exec(Option_t *)
        ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
        ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
        
-       if(ch1 == ch2){
+       if(ch1 == ch2 && !fGeneratorOnly){
          if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
            fPairSplitCut[3][i]->AddAt('1',j);
            continue;
@@ -3165,6 +3212,7 @@ void AliChaoticity::Exec(Option_t *)
                MomResCorr12 = fMomResC2->GetBinContent(momResIndex+1, momBin12);
                MomResCorr13 = fMomResC2->GetBinContent(momResIndex+1, momBin13);
                MomResCorr23 = fMomResC2->GetBinContent(momResIndex+1, momBin23);
+               
                if(MomResCorr12 > 1.2 || MomResCorr13 > 1.2 || MomResCorr23 > 1.2) {// Safety check
                  if(fMbin==0 && bin1==0) {
                    ((TH3F*)fOutputList->FindObject("fTPNRejects3"))->Fill(qinv12, qinv13, qinv23, sqrt(fabs(weight12*weight13*weight23)));
@@ -3865,7 +3913,7 @@ void AliChaoticity::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt,
   qSide = fabs(qSide);
   qLong = fabs(qLong);
   Float_t wd=0, xd=0, yd=0, zd=0;
-  //Float_t qinv_temp=GetQinv(0,track1, track2);
+  //Float_t qinvtemp=GetQinv(0,track1, track2);
   //
   
   if(kt < fKmeanT[0]) {fKtIndexL=0; fKtIndexH=1;}
index 4fb4b70..95eded7 100644 (file)
@@ -81,6 +81,7 @@ class AliChaoticity : public AliAnalysisTaskSE {
   void SetTabulatePairs(Bool_t tabulate) {fTabulatePairs = tabulate;}
   void SetPbPbCase(Bool_t pbpb) {fPbPbcase = pbpb;}
   void SetGenerateSignal(Bool_t gen) {fGenerateSignal = gen;}
+  void SetGeneratorOnly(Bool_t genOnly) {fGeneratorOnly = genOnly;}
   void SetCentBinRange(Int_t low, Int_t high) {fCentBinLowLimit = low; fCentBinHighLimit = high;}
   void SetLEGOCase(Bool_t lego) {fLEGO = lego;}
   void SetFilterBit(UInt_t filterbit) {fFilterBit = filterbit;}
@@ -235,6 +236,7 @@ class AliChaoticity : public AliAnalysisTaskSE {
   Bool_t fAODcase;
   Bool_t fPbPbcase;
   Bool_t fGenerateSignal;
+  Bool_t fGeneratorOnly;
   Bool_t fPdensityExplicitLoop;
   Bool_t fPdensityPairCut;
   Bool_t fTabulatePairs;