In AOD Event Reader: TOF mismatch with new method + macros for proton femtoscopy...
authormajanik <majanik@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 17 Jun 2013 15:33:13 +0000 (15:33 +0000)
committermajanik <majanik@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 17 Jun 2013 15:33:13 +0000 (15:33 +0000)
PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx
PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.cxx
PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.h
PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/AvgSep/ConfigFemtoAnalysis.C [new file with mode: 0644]
PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcOnly/ConfigFemtoAnalysis.C [new file with mode: 0644]
PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcOnly/monitors/ConfigFemtoAnalysis.C [new file with mode: 0644]
PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcTof/central/ConfigFemtoAnalysis.C [new file with mode: 0644]
PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcTof/monitors/ConfigFemtoAnalysis.C [new file with mode: 0644]
PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcTof/semicentral/ConfigFemtoAnalysis.C [new file with mode: 0644]

index 755f1ec21ab5636fc24c6ad4cdd4906a5b049c02..538bca8096b0a93c8a50c81cf5d341ef3e16f18c 100644 (file)
@@ -32,7 +32,7 @@
 ClassImp(AliFemtoEventReaderAOD)
 
 #if !(ST_NO_NAMESPACES)
-  using namespace units;
+using namespace units;
 #endif
 
 using namespace std;
@@ -40,7 +40,7 @@ using namespace std;
 double fV1[3];
 
 //____________________________
-//constructor with 0 parameters , look at default settings 
+//constructor with 0 parameters , look at default settings
 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
   fNumberofEvent(0),
   fCurEvent(0),
@@ -132,7 +132,7 @@ AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventRea
 {
   // assignment operator
   if (this == &aReader)
-   return *this;
+    return *this;
 
   fInputFile = aReader.fInputFile;
   fNumberofEvent = aReader.fNumberofEvent;
@@ -168,7 +168,7 @@ AliFemtoString AliFemtoEventReaderAOD::Report()
 //__________________
 void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
 {
-  //setting the name of file where names of AOD file are written 
+  //setting the name of file where names of AOD file are written
   //it takes only this files which have good trees
   char buffer[256];
   fInputFile=string(inputFile);
@@ -177,26 +177,26 @@ void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
   fTree = new TChain("aodTree");
 
   if(infile.good()==true)
-    { 
-      //checking if all give files have good tree inside
-      while (infile.eof()==false)
-       {
-         infile.getline(buffer,256);
-         TFile *aodFile=TFile::Open(buffer,"READ");
-         if (aodFile!=0x0)
-           {   
+  {
+    //checking if all give files have good tree inside
+    while (infile.eof()==false)
+    {
+      infile.getline(buffer,256);
+      TFile *aodFile=TFile::Open(buffer,"READ");
+      if (aodFile!=0x0)
+           {
              TTree* tree = (TTree*) aodFile->Get("aodTree");
              if (tree!=0x0)
-               {
-                 //              cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
-                 fTree->AddFile(buffer);
-                 delete tree;
-               }
-             aodFile->Close(); 
+        {
+          //             cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
+          fTree->AddFile(buffer);
+          delete tree;
+        }
+             aodFile->Close();
            }
-         delete aodFile;
-       }
+      delete aodFile;
     }
+  }
 }
 
 AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
@@ -205,15 +205,15 @@ AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
   // convert it to AliFemtoEvent and return
   // for further analysis
   AliFemtoEvent *hbtEvent = 0;
-    // cout<<"reader"<<endl;
-  if (fCurEvent==fNumberofEvent)//open next file  
+  // cout<<"reader"<<endl;
+  if (fCurEvent==fNumberofEvent)//open next file
+  {
+    if(fNumberofEvent==0)
     {
-      if(fNumberofEvent==0)    
-       {
-         fEvent=new AliAODEvent();
-         fEvent->ReadFromTree(fTree);
+      fEvent=new AliAODEvent();
+      fEvent->ReadFromTree(fTree);
 
-         // Check for the existence of the additional information
+      // Check for the existence of the additional information
 //       fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
 
 //       if (fPWG2AODTracks) {
@@ -221,29 +221,29 @@ AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
 //         cout << "Reading only tracks with the additional information" << endl;
 //       }
 
-         fNumberofEvent=fTree->GetEntries();
-         //      cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
-         fCurEvent=0;
-       }
-      else //no more data to read
-       {
-            // cout<<"no more files "<<hbtEvent<<endl;
-         fReaderStatus=1;
-         return hbtEvent; 
-       }
-    }          
+      fNumberofEvent=fTree->GetEntries();
+      //         cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
+      fCurEvent=0;
+    }
+    else //no more data to read
+    {
+      // cout<<"no more files "<<hbtEvent<<endl;
+      fReaderStatus=1;
+      return hbtEvent;
+    }
+  }
 
-    // cout<<"starting to read event "<<fCurEvent<<endl;
+  // cout<<"starting to read event "<<fCurEvent<<endl;
   fTree->GetEvent(fCurEvent);//getting next event
   //  cout << "Read event " << fEvent << " from file " << fTree << endl;
-       
+
   hbtEvent = new AliFemtoEvent;
 
   CopyAODtoFemtoEvent(hbtEvent);
   fCurEvent++;
 
 
-  return hbtEvent; 
+  return hbtEvent;
 }
 
 void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
@@ -264,7 +264,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
   tEvent->SetZDCParticipants(0);
   tEvent->SetTriggerMask(fEvent->GetTriggerMask());
   tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
-  
+
   // Attempt to access MC header
   AliAODMCHeader *mcH;
   TClonesArray *mcP=0;
@@ -292,37 +292,37 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
     for (int ip=0; ip<mcP->GetEntries(); ip++) {
       motherpart = (AliAODMCParticle *) mcP->At(ip);
       if (motherpart->GetDaughter(0) > 0)
-       motherids[motherpart->GetDaughter(0)] = ip;
+        motherids[motherpart->GetDaughter(0)] = ip;
       if (motherpart->GetDaughter(1) > 0)
-       motherids[motherpart->GetDaughter(1)] = ip;
+        motherids[motherpart->GetDaughter(1)] = ip;
     }
   }
 
   // Primary Vertex position
   //  double fV1[3];
   if(fpA2013)
-    {
-      const AliAODVertex* trkVtx = fEvent->GetPrimaryVertex();
-      if (!trkVtx || trkVtx->GetNContributors()<=0)  return;
-      TString vtxTtl = trkVtx->GetTitle();
-      if (!vtxTtl.Contains("VertexerTracks")) return;
-      Float_t zvtx = trkVtx->GetZ();
-      const AliAODVertex* spdVtx = fEvent->GetPrimaryVertexSPD();
-      if (spdVtx->GetNContributors()<=0)  return;
-      TString vtxTyp = spdVtx->GetTitle();
-      Double_t cov[6]={0};
-      spdVtx->GetCovarianceMatrix(cov);
-      Double_t zRes = TMath::Sqrt(cov[5]);
-      if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25))  return;
-      if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)  return;
-
-      if (TMath::Abs(zvtx) > 10)  return;
-    }
+  {
+    const AliAODVertex* trkVtx = fEvent->GetPrimaryVertex();
+    if (!trkVtx || trkVtx->GetNContributors()<=0)  return;
+    TString vtxTtl = trkVtx->GetTitle();
+    if (!vtxTtl.Contains("VertexerTracks")) return;
+    Float_t zvtx = trkVtx->GetZ();
+    const AliAODVertex* spdVtx = fEvent->GetPrimaryVertexSPD();
+    if (spdVtx->GetNContributors()<=0)  return;
+    TString vtxTyp = spdVtx->GetTitle();
+    Double_t cov[6]={0};
+    spdVtx->GetCovarianceMatrix(cov);
+    Double_t zRes = TMath::Sqrt(cov[5]);
+    if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25))  return;
+    if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)  return;
+
+    if (TMath::Abs(zvtx) > 10)  return;
+  }
   fEvent->GetPrimaryVertex()->GetPosition(fV1);
 
   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
   tEvent->SetPrimVertPos(vertex);
-       
+
   //starting to reading tracks
   int nofTracks=0;  //number of reconstructed tracks in event
 
@@ -331,32 +331,32 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
   //     nofTracks=fPWG2AODTracks->GetEntries();
   //   else
   nofTracks=fEvent->GetNumberOfTracks();
-  
+
   AliEventplane *ep = fEvent->GetEventplane();
   if (ep) {
     tEvent->SetEP(ep);
-        if (fisEPVZ)
-            tEvent->SetReactionPlaneAngle(ep->GetEventplane("V0",fEvent,2));
-        else
-    tEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
+    if (fisEPVZ)
+      tEvent->SetReactionPlaneAngle(ep->GetEventplane("V0",fEvent,2));
+    else
+      tEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
   }
 
   AliCentrality *cent = fEvent->GetCentrality();
-  
+
   if (!fEstEventMult && cent && fUsePreCent) {
     if ((cent->GetCentralityPercentile("V0M")*10 < fCentRange[0]) ||
-       (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
-      {
-            // cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
-       
-       return;
-      }
+        (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
+    {
+      // cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
+
+      return;
+    }
   }
 
   int realnofTracks=0;   // number of track which we use in a analysis
-  int tracksPrim=0;     
+  int tracksPrim=0;
 
-  int labels[20000];   
+  int labels[20000];
   for (int il=0; il<20000; il++) labels[il] = -1;
 
   // looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
@@ -370,46 +370,46 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 
   int tNormMult = 0;
   for (int i=0;i<nofTracks;i++)
-    {
-      AliFemtoTrack* trackCopy = new AliFemtoTrack();  
+  {
+    AliFemtoTrack* trackCopy = new AliFemtoTrack();
 
 //       if (fPWG2AODTracks) {
 //     // Read tracks from the additional pwg2 specific AOD part
 //     // if they exist
-//     // Note that in that case all the AOD tracks without the 
+//     // Note that in that case all the AOD tracks without the
 //     // additional information will be ignored !
 //     AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
 
 //     // Getting the AOD track through the ref of the additional info
-//     AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack(); 
+//     AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
 //     if (!aodtrack->TestFilterBit(fFilterBit)) {
 //       delete trackCopy;
 //       continue;
 //     }
 
+
 //     if (aodtrack->IsOn(AliESDtrack::kTPCrefit))
-//       if (aodtrack->Chi2perNDF() < 6.0) 
+//       if (aodtrack->Chi2perNDF() < 6.0)
 //         if (aodtrack->Eta() < 0.9)
 //           tNormMult++;
 
 
 //     CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
-       
+
 //     if (mcP) {
 //       // Fill the hidden information with the simulated data
 //       //      Int_t pLabel = aodtrack->GetLabel();
 //       AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
 
 //       // Check the mother information
-         
+
 //       // Using the new way of storing the freeze-out information
 //       // Final state particle is stored twice on the stack
 //       // one copy (mother) is stored with original freeze-out information
 //       //   and is not tracked
 //       // the other one (daughter) is stored with primary vertex position
 //       //   and is tracked
-         
+
 //       // Freeze-out coordinates
 //       double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
 //       fpx = tPart->Xv() - fV1[0];
@@ -424,7 +424,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 //       fpy *= 1e13;
 //       fpz *= 1e13;
 //       fpt *= 1e13;
-         
+
 //       //      cout << "Looking for mother ids " << endl;
 //       if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
 //         //  cout << "Got mother id" << endl;
@@ -434,27 +434,27 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 //           // It is the same particle
 //           // Read in the original freeze-out information
 //           // and convert it from to [fm]
-             
-//           // EPOS style 
+
+//           // EPOS style
 //           //          fpx = mother->Xv()*1e13*0.197327;
 //           //          fpy = mother->Yv()*1e13*0.197327;
 //           //          fpz = mother->Zv()*1e13*0.197327;
 //           //          fpt = mother->T() *1e13*0.197327*0.5;
-             
-             
-//           // Therminator style 
+
+
+//           // Therminator style
 //           fpx = mother->Xv()*1e13;
 //           fpy = mother->Yv()*1e13;
 //           fpz = mother->Zv()*1e13;
 //           fpt = mother->T() *1e13*3e10;
-             
+
 //         }
 //       }
-         
+
 //       //       if (fRotateToEventPlane) {
 //       //    double tPhi = TMath::ATan2(fpy, fpx);
 //       //    double tRad = TMath::Hypot(fpx, fpy);
-       
+
 //       //    fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
 //       //    fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
 //       //       }
@@ -464,7 +464,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 //       //      if (fRotateToEventPlane) {
 //       //        double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
 //       //        double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
-           
+
 //       //        tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
 //       //                               tRad*TMath::Sin(tPhi - tReactionPlane),
 //       //                               tPart->Pz());
@@ -477,9 +477,9 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 //                         tPart->Pz()*tPart->Pz());
 //       if (mass2>0.0)
 //         tInfo->SetMass(TMath::Sqrt(mass2));
-//       else 
+//       else
 //         tInfo->SetMass(0.0);
-         
+
 //       tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
 //       trackCopy->SetHiddenInfo(tInfo);
 
@@ -495,174 +495,183 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 //     }
 //       }
 //       else {
-       // No additional information exists
-       // Read in the normal AliAODTracks 
+    // No additional information exists
+    // Read in the normal AliAODTracks
 
-       //      const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
-       AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
-       
+    // const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
+    AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
 
 
-       if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
-       
-       if (fFilterBit && !aodtrack->TestFilterBit(fFilterBit)) {
-         delete trackCopy;
-         continue;
-       }
 
-       if (fFilterMask && !aodtrack->TestFilterBit(fFilterMask)) {
-         delete trackCopy;
-         continue;
-       }               
+    if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
 
-       //counting particles to set multiplicity
-       double impact[2];
-       double covimpact[3];
-       if (aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
-         if(impact[0]<0.2 && TMath::Abs(impact[1]+fV1[2])<2.0)
-           //if (aodtrack->IsPrimaryCandidate()) //? instead of kinks?
-             if (aodtrack->Chi2perNDF() < 4.0) 
-               if (aodtrack->Pt() > 0.15 && aodtrack->Pt() < 20) 
-                 if (aodtrack->GetTPCNcls() > 70)
-                   if (aodtrack->Eta() < 0.8)
-                     tNormMult++;
-       } 
-
-       CopyAODtoFemtoTrack(aodtrack, trackCopy);
-
-       // copying PID information from the correspondent track
-       //      const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
-
-
-       AliAODTrack *aodtrackpid;
-       if((fFilterBit ==  (1 << (7))) || fFilterMask==128) //for TPC Only tracks we have to copy PID information from corresponding global tracks
-         aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
-       else
-         aodtrackpid = fEvent->GetTrack(i);
-        CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
-       
-       if (mcP) {
-         // Fill the hidden information with the simulated data
-         //      Int_t pLabel = aodtrack->GetLabel();
-         AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
-         
-         AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
-         double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
-         if (!tPart) {
-           fpx = fV1[0];
-           fpy = fV1[1];
-           fpz = fV1[2];
-           tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
-           tInfo->SetPDGPid(0);
-           tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
-           tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
-           tInfo->SetMass(0);
-         }
-         else {
-           // Check the mother information
-         
-           // Using the new way of storing the freeze-out information
-           // Final state particle is stored twice on the stack
-           // one copy (mother) is stored with original freeze-out information
-           //   and is not tracked
-           // the other one (daughter) is stored with primary vertex position
-           //   and is tracked
-           
-           // Freeze-out coordinates
-           fpx = tPart->Xv() - fV1[0];
-           fpy = tPart->Yv() - fV1[1];
-           fpz = tPart->Zv() - fV1[2];
-           //    fpt = tPart->T();
-           
-           tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
-           
-           fpx *= 1e13;
-           fpy *= 1e13;
-           fpz *= 1e13;
-           //    fpt *= 1e13;
-           
-           //      cout << "Looking for mother ids " << endl;
-           if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
-             //        cout << "Got mother id" << endl;
-             AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
-             // Check if this is the same particle stored twice on the stack
-             if (mother) {
-               if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
-                 // It is the same particle
-                 // Read in the original freeze-out information
-                 // and convert it from to [fm]
-                 
-                 // EPOS style 
-                 //      fpx = mother->Xv()*1e13*0.197327;
-                 //      fpy = mother->Yv()*1e13*0.197327;
-                 //      fpz = mother->Zv()*1e13*0.197327;
-                 //      fpt = mother->T() *1e13*0.197327*0.5;
-                 
-                 
-                 // Therminator style 
-                 fpx = mother->Xv()*1e13;
-                 fpy = mother->Yv()*1e13;
-                 fpz = mother->Zv()*1e13;
-                 //          fpt = mother->T() *1e13*3e10;
-                 
-               }
-             }
-           }
-           
-           //       if (fRotateToEventPlane) {
-           //  double tPhi = TMath::ATan2(fpy, fpx);
-           //  double tRad = TMath::Hypot(fpx, fpy);
-           
-           //  fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
-           //  fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
-           //       }
-           
-           tInfo->SetPDGPid(tPart->GetPdgCode());
-           
-           //    if (fRotateToEventPlane) {
-           //      double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
-           //      double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
-           
-           //      tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
-           //                             tRad*TMath::Sin(tPhi - tReactionPlane),
-           //                             tPart->Pz());
-           //    }
-           //       else
-           tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
-           Double_t mass2 = (tPart->E() *tPart->E() -
-                             tPart->Px()*tPart->Px() -
-                             tPart->Py()*tPart->Py() -
-                             tPart->Pz()*tPart->Pz());
-           if (mass2>0.0)
-             tInfo->SetMass(TMath::Sqrt(mass2));
-           else 
-             tInfo->SetMass(0.0);
-           
-           tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
-         }
-         trackCopy->SetHiddenInfo(tInfo);
-       }
+    if (fFilterBit && !aodtrack->TestFilterBit(fFilterBit)) {
+      delete trackCopy;
+      continue;
+    }
 
-       double pxyz[3];
+    if (fFilterMask && !aodtrack->TestFilterBit(fFilterMask)) {
+      delete trackCopy;
+      continue;
+    }
 
-       //AliExternalTrackParam *param = new AliExternalTrackParam(*aodtrack->GetInnerParam());
-       trackCopy->SetInnerMomentum(aodtrack->GetTPCmomentum());
+    // cout << "Muon? " << aodtrack->IsMuonTrack() << endl;
+    // cout << "Type? " << aodtrack->GetType() << endl;
+
+    // if (aodtrack->IsMuonTrack()) {
+    //   cout << "muon" << endl;
+    //   delete trackCopy;
+    //   continue;
+    // }
+
+    //counting particles to set multiplicity
+    double impact[2];
+    double covimpact[3];
+    if (aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
+      if(impact[0]<0.2 && TMath::Abs(impact[1]+fV1[2])<2.0)
+        //if (aodtrack->IsPrimaryCandidate()) //? instead of kinks?
+             if (aodtrack->Chi2perNDF() < 4.0)
+          if (aodtrack->Pt() > 0.15 && aodtrack->Pt() < 20)
+            if (aodtrack->GetTPCNcls() > 70)
+              if (aodtrack->Eta() < 0.8)
+                tNormMult++;
+    }
 
-       aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
-       const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
-       // Check the sanity of the tracks - reject zero momentum tracks
-       if (ktP.Mag() == 0) {
-         delete trackCopy;
-         continue;
-       }
-       //    }
-  
-       
-       tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
-       realnofTracks++;//real number of tracks         
+    CopyAODtoFemtoTrack(aodtrack, trackCopy);
+
+    // copying PID information from the correspondent track
+    // const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
+
+
+    AliAODTrack *aodtrackpid;
+    if((fFilterBit ==  (1 << (7))) || fFilterMask==128) //for TPC Only tracks we have to copy PID information from corresponding global tracks
+      aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
+    else
+      aodtrackpid = fEvent->GetTrack(i);
+    CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
+
+    if (mcP) {
+      // Fill the hidden information with the simulated data
+      //         Int_t pLabel = aodtrack->GetLabel();
+      AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
+
+      AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
+      double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
+      if (!tPart) {
+        fpx = fV1[0];
+        fpy = fV1[1];
+        fpz = fV1[2];
+        tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
+        tInfo->SetPDGPid(0);
+        tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
+        tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
+        tInfo->SetMass(0);
+      }
+      else {
+        // Check the mother information
+
+        // Using the new way of storing the freeze-out information
+        // Final state particle is stored twice on the stack
+        // one copy (mother) is stored with original freeze-out information
+        //   and is not tracked
+        // the other one (daughter) is stored with primary vertex position
+        //   and is tracked
+
+        // Freeze-out coordinates
+        fpx = tPart->Xv() - fV1[0];
+        fpy = tPart->Yv() - fV1[1];
+        fpz = tPart->Zv() - fV1[2];
+        //       fpt = tPart->T();
+
+        tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
+
+        fpx *= 1e13;
+        fpy *= 1e13;
+        fpz *= 1e13;
+        //       fpt *= 1e13;
+
+        //      cout << "Looking for mother ids " << endl;
+        if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
+          //   cout << "Got mother id" << endl;
+          AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
+          // Check if this is the same particle stored twice on the stack
+          if (mother) {
+            if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
+              // It is the same particle
+              // Read in the original freeze-out information
+              // and convert it from to [fm]
+
+              // EPOS style
+              //         fpx = mother->Xv()*1e13*0.197327;
+              //         fpy = mother->Yv()*1e13*0.197327;
+              //         fpz = mother->Zv()*1e13*0.197327;
+              //         fpt = mother->T() *1e13*0.197327*0.5;
+
+
+              // Therminator style
+              fpx = mother->Xv()*1e13;
+              fpy = mother->Yv()*1e13;
+              fpz = mother->Zv()*1e13;
+              //             fpt = mother->T() *1e13*3e10;
+
+            }
+          }
+        }
+
+        //       if (fRotateToEventPlane) {
+        //     double tPhi = TMath::ATan2(fpy, fpx);
+        //     double tRad = TMath::Hypot(fpx, fpy);
+
+        //     fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
+        //     fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
+        //       }
+
+        tInfo->SetPDGPid(tPart->GetPdgCode());
+
+        //       if (fRotateToEventPlane) {
+        //         double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
+        //         double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
+
+        //         tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
+        //                                tRad*TMath::Sin(tPhi - tReactionPlane),
+        //                                tPart->Pz());
+        //       }
+        //       else
+        tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
+        Double_t mass2 = (tPart->E() *tPart->E() -
+                          tPart->Px()*tPart->Px() -
+                          tPart->Py()*tPart->Py() -
+                          tPart->Pz()*tPart->Pz());
+        if (mass2>0.0)
+          tInfo->SetMass(TMath::Sqrt(mass2));
+        else
+          tInfo->SetMass(0.0);
+
+        tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
+      }
+      trackCopy->SetHiddenInfo(tInfo);
     }
-  
-  tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event    
+
+    double pxyz[3];
+
+    //AliExternalTrackParam *param = new AliExternalTrackParam(*aodtrack->GetInnerParam());
+    trackCopy->SetInnerMomentum(aodtrack->GetTPCmomentum());
+
+    aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
+    const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
+    // Check the sanity of the tracks - reject zero momentum tracks
+    if (ktP.Mag() == 0) {
+      delete trackCopy;
+      continue;
+    }
+    //    }
+
+
+    tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
+    realnofTracks++;//real number of tracks
+  }
+
+  tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
   tEvent->SetNormalizedMult(tracksPrim);
 
   if (cent) {
@@ -696,7 +705,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
   else if (fEstEventMult==kCentralityZNA) {
     if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNA")));
   }
- else if (fEstEventMult==kCentralityZNC) {
 else if (fEstEventMult==kCentralityZNC) {
     if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNC")));
   }
   else if (fEstEventMult==kCentralityCL1) {
@@ -717,90 +726,90 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
   else if (fEstEventMult==kCentralityNPA) {
     if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("NPA")));
   }
- else if (fEstEventMult==kCentralityFMD) {
 else if (fEstEventMult==kCentralityFMD) {
     if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("FMD")));
   }
   else if(fEstEventMult==kGlobalCount){
     tEvent->SetNormalizedMult(tNormMult); //particles counted in the loop, trying to reproduce GetReferenceMultiplicity. If better (default) method appears it should be changed
   }
   else if(fEstEventMult==kReference)
-    {
-      tEvent->SetNormalizedMult(fAODheader->GetRefMultiplicity());
-    }
+  {
+    tEvent->SetNormalizedMult(fAODheader->GetRefMultiplicity());
+  }
   else if(fEstEventMult==kTPCOnlyRef)
-    {
-      tEvent->SetNormalizedMult(fAODheader->GetTPConlyRefMultiplicity());
-    }
+  {
+    tEvent->SetNormalizedMult(fAODheader->GetTPConlyRefMultiplicity());
+  }
   else if(fEstEventMult == kVZERO)
-    {
-      Float_t multV0 = 0;
-      for (Int_t i=0; i<64; i++)
-       multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
-      tEvent->SetNormalizedMult(multV0);
-    }
+  {
+    Float_t multV0 = 0;
+    for (Int_t i=0; i<64; i++)
+      multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
+    tEvent->SetNormalizedMult(multV0);
+  }
 
   if (mcP) delete [] motherids;
 
-    // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
+  // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 
   if(fReadV0)
-    {
-      int count_pass = 0;
-      for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
-       AliAODv0* aodv0 = fEvent->GetV0(i);
-       if (!aodv0) continue;
-       if(aodv0->GetNDaughters()>2) continue;
-       if(aodv0->GetNProngs()>2) continue;
-       if(aodv0->GetCharge()!=0) continue;
-       if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
-       if(aodv0->CosPointingAngle(fV1)<0.998) continue;
-       AliFemtoV0* trackCopyV0 = new AliFemtoV0();
-       count_pass++;
-       CopyAODtoFemtoV0(aodv0, trackCopyV0);
-       tEvent->V0Collection()->push_back(trackCopyV0);
-       //cout<<"Pushback v0 to v0collection"<<endl;
-      }
+  {
+    int count_pass = 0;
+    for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
+      AliAODv0* aodv0 = fEvent->GetV0(i);
+      if (!aodv0) continue;
+      if(aodv0->GetNDaughters()>2) continue;
+      if(aodv0->GetNProngs()>2) continue;
+      if(aodv0->GetCharge()!=0) continue;
+      if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
+      if(aodv0->CosPointingAngle(fV1)<0.998) continue;
+      AliFemtoV0* trackCopyV0 = new AliFemtoV0();
+      count_pass++;
+      CopyAODtoFemtoV0(aodv0, trackCopyV0);
+      tEvent->V0Collection()->push_back(trackCopyV0);
+      //cout<<"Pushback v0 to v0collection"<<endl;
     }
+  }
 
 }
 
-void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack, 
-                                                AliFemtoTrack *tFemtoTrack 
-                                                //                                              AliPWG2AODTrack *tPWG2AODTrack
-                                                )
+void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
+                                                 AliFemtoTrack *tFemtoTrack
+                                                 //                                             AliPWG2AODTrack *tPWG2AODTrack
+  )
 {
   // Copy the track information from the AOD into the internal AliFemtoTrack
   // If it exists, use the additional information from the PWG2 AOD
 
   // Primary Vertex position
-  
+
   fEvent->GetPrimaryVertex()->GetPosition(fV1);
   //  fEvent->GetPrimaryVertex()->GetXYZ(fV1);
 
   tFemtoTrack->SetCharge(tAodTrack->Charge());
-  
+
   double pxyz[3];
   tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
   AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
   tFemtoTrack->SetP(v);//setting momentum
   tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
   const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
-  //setting track helix 
+  //setting track helix
   const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
-  AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge())); 
+  AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
   tFemtoTrack->SetHelix(helix);
-               
+
   // Flags
   tFemtoTrack->SetTrackId(tAodTrack->GetID());
   tFemtoTrack->SetFlags(tAodTrack->GetFlags());
   tFemtoTrack->SetLabel(tAodTrack->GetLabel());
-               
-  // Track quality information 
+
+  // Track quality information
   float covmat[6];
-  tAodTrack->GetCovMatrix(covmat);  
+  tAodTrack->GetCovMatrix(covmat);
+
+  // ! DCA information is done in CopyPIDtoFemtoTrack()
 
-        // ! DCA information is done in CopyPIDtoFemtoTrack()
-  
        // double impact[2];
        // double covimpact[3];
 
@@ -820,7 +829,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
   //   else
   //     tFemtoTrack->SetImpactD(0.0);
   //   tFemtoTrack->SetImpactD(tAodTrack->DCA());
-    
+
   //   tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
 
 
@@ -828,18 +837,18 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
   //   tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
 
 
-  //   cout 
-    //    << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]) 
-    //    << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0] 
-    //    << tAodTrack->Yv() - fV1[1] 
-//     << "xv = " << tAodTrack->Xv() << endl 
-//     << "fv1[0] = " << fV1[0]  << endl 
-//     << "yv = " << tAodTrack->Yv()  << endl 
-//     << "fv1[1] = " << fV1[1]  << endl 
-//     << "zv = " << tAodTrack->Zv()  << endl 
-//     << "fv1[2] = " << fV1[2]  << endl 
-//     << "impact[0] = " << impact[0]  << endl 
-//     << "impact[1] = " << impact[1]  << endl 
+  //   cout
+  //    << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
+  //    << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
+  //    << tAodTrack->Yv() - fV1[1]
+//     << "xv = " << tAodTrack->Xv() << endl
+//     << "fv1[0] = " << fV1[0]  << endl
+//     << "yv = " << tAodTrack->Yv()  << endl
+//     << "fv1[1] = " << fV1[1]  << endl
+//     << "zv = " << tAodTrack->Zv()  << endl
+//     << "fv1[2] = " << fV1[2]  << endl
+//     << "impact[0] = " << impact[0]  << endl
+//     << "impact[1] = " << impact[1]  << endl
 //     << endl << endl ;
 
   tFemtoTrack->SetCdd(covmat[0]);
@@ -850,15 +859,15 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
   tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
   tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
-  tFemtoTrack->SetTPCsignalN(1); 
-  tFemtoTrack->SetTPCsignalS(1); 
+  tFemtoTrack->SetTPCsignalN(1);
+  tFemtoTrack->SetTPCsignalS(1);
   tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
 
 //   if (tPWG2AODTrack) {
 //     // Copy the PWG2 specific information if it exists
 //     tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
 //     tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
-    
+
 //     double xtpc[3] = {0,0,0};
 //     tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
 //     tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
@@ -866,10 +875,10 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
 //     tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 //   }
 //   else {
-    // If not use dummy values
+  // If not use dummy values
   tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
   tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
-  
+
 
   float globalPositionsAtRadii[9][3];
   float bfield = 5*fMagFieldSign;
@@ -881,22 +890,22 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
     tpcPositions[i] = new double[3];
   double tpcExit[3]={globalPositionsAtRadii[8][0],globalPositionsAtRadii[8][1],globalPositionsAtRadii[8][2]};
   for(int i=0;i<9;i++)
-    {
-      tpcPositions[i][0] = globalPositionsAtRadii[i][0];
-      tpcPositions[i][1] = globalPositionsAtRadii[i][1];
-      tpcPositions[i][2] = globalPositionsAtRadii[i][2];
-    }
+  {
+    tpcPositions[i][0] = globalPositionsAtRadii[i][0];
+    tpcPositions[i][1] = globalPositionsAtRadii[i][1];
+    tpcPositions[i][2] = globalPositionsAtRadii[i][2];
+  }
   tFemtoTrack->SetNominalTPCEntrancePoint(tpcEntrance);
   tFemtoTrack->SetNominalTPCPoints(tpcPositions);
   tFemtoTrack->SetNominalTPCExitPoint(tpcExit);
-    for(int i=0;i<9;i++)
-      delete [] tpcPositions[i];
-    delete [] tpcPositions;
+  for(int i=0;i<9;i++)
+    delete [] tpcPositions[i];
+  delete [] tpcPositions;
   //   }
-  
+
   //   //  cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << "  " << tAodTrack->Zv() << "  " << tAodTrack->GetTPCNcls() << endl;
-  
-  
+
+
   int indexes[3];
   for (int ik=0; ik<3; ik++) {
     indexes[ik] = 0;
@@ -935,8 +944,8 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
   tFemtoV0->SetmomNeg(momneg);
 
   //jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
-  //void SettpcHitsPos(const int& i);      
-  //void SettpcHitsNeg(const int& i);      
+  //void SettpcHitsPos(const int& i);
+  //void SettpcHitsNeg(const int& i);
 
   //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
   //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
@@ -957,10 +966,10 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
   tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
   tFemtoV0->SetrapLambda(tAODv0->RapLambda());
   tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
-  
-  //void SetcTauLambda( float x);   
-  //void SetcTauK0Short( float x); 
-  
+
+  //void SetcTauLambda( float x);
+  //void SetcTauK0Short( float x);
+
   //tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
   tFemtoV0->SetptV0(tAODv0->Pt());
   tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
@@ -968,7 +977,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
   //tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
   //tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
   //tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
-  
+
   tFemtoV0->SetidNeg(tAODv0->GetNegID());
   //cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
   //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
@@ -993,78 +1002,78 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
   AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
 
   if(trackpos && trackneg)
+  {
+    tFemtoV0->SetEtaPos(trackpos->Eta());
+    tFemtoV0->SetEtaNeg(trackneg->Eta());
+    tFemtoV0->SetptotPos(tAODv0->PProng(0));
+    tFemtoV0->SetptotNeg(tAODv0->PProng(1));
+    tFemtoV0->SetptPos(trackpos->Pt());
+    tFemtoV0->SetptNeg(trackneg->Pt());
+
+    //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
+    //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
+    tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
+    tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
+    tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
+    tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
+    tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
+    tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
+    tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
+    tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
+    tFemtoV0->SetStatusPos(trackpos->GetStatus());
+    tFemtoV0->SetStatusNeg(trackneg->GetStatus());
+
+    tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
+    tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
+    tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
+    tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
+    tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
+    tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
+
+
+    float bfield = 5*fMagFieldSign;
+    float globalPositionsAtRadiiPos[9][3];
+    GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
+    double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
+    double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
+
+    float globalPositionsAtRadiiNeg[9][3];
+    GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
+    double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
+    double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
+
+    AliFemtoThreeVector tmpVec;
+    tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
+    tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
+
+    tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
+    tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
+
+    tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
+    tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
+
+    tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
+    tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
+
+    AliFemtoThreeVector vecTpcPos[9];
+    AliFemtoThreeVector vecTpcNeg[9];
+    for(int i=0;i<9;i++)
     {
-      tFemtoV0->SetEtaPos(trackpos->Eta());
-      tFemtoV0->SetEtaNeg(trackneg->Eta());
-      tFemtoV0->SetptotPos(tAODv0->PProng(0));
-      tFemtoV0->SetptotNeg(tAODv0->PProng(1));
-      tFemtoV0->SetptPos(trackpos->Pt());
-      tFemtoV0->SetptNeg(trackneg->Pt());
-
-      //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
-      //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
-      tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
-      tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
-      tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
-      tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
-      tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
-      tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
-      tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
-      tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
-      tFemtoV0->SetStatusPos(trackpos->GetStatus());
-      tFemtoV0->SetStatusNeg(trackneg->GetStatus());
-
-      tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
-      tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
-      tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
-      tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
-      tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
-      tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
-
-
-      float bfield = 5*fMagFieldSign;
-      float globalPositionsAtRadiiPos[9][3];
-      GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
-      double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
-      double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
-
-      float globalPositionsAtRadiiNeg[9][3];
-      GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
-      double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
-      double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
-
-      AliFemtoThreeVector tmpVec;
-      tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
-      tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
-
-      tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
-      tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
-
-      tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
-      tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
-
-      tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
-      tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
-
-      AliFemtoThreeVector vecTpcPos[9];
-      AliFemtoThreeVector vecTpcNeg[9];
-      for(int i=0;i<9;i++)
-       {
-         vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
-         vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
-       }
-      tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
-      tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
+      vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
+      vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
+    }
+    tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
+    tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
 
-      tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
-      tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
+    tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
+    tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
 
-      tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
-      tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
+    tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
+    tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
 
-        if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 )
-       {
-            if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 )
+    if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 )
+    {
+      if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 )
            {
              tFemtoV0->SetPosNSigmaTOFK(-1000);
              tFemtoV0->SetNegNSigmaTOFK(-1000);
@@ -1080,36 +1089,36 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFem
              tFemtoV0->SetTOFPionTimeNeg(-1000);
              tFemtoV0->SetTOFKaonTimeNeg(-1000);
            }
-       }
-      else
-       {
-         tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
-         tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
-         tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
-         tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
-         tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
-         tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
-
-         double TOFSignalPos = trackpos->GetTOFsignal();
-         double TOFSignalNeg = trackneg->GetTOFsignal();
-         double pidPos[5];
-         double pidNeg[5];
-         trackpos->GetIntegratedTimes(pidPos);
-         trackneg->GetIntegratedTimes(pidNeg);
-
-         tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
-         tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
-         tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
-         tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
-         tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
-         tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
-       }
     }
-  else
+    else
     {
-      tFemtoV0->SetStatusPos(999);
-      tFemtoV0->SetStatusNeg(999);
+      tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
+      tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
+      tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
+      tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
+      tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
+      tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
+
+      double TOFSignalPos = trackpos->GetTOFsignal();
+      double TOFSignalNeg = trackneg->GetTOFsignal();
+      double pidPos[5];
+      double pidNeg[5];
+      trackpos->GetIntegratedTimes(pidPos);
+      trackneg->GetIntegratedTimes(pidNeg);
+
+      tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
+      tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
+      tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
+      tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
+      tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
+      tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
     }
+  }
+  else
+  {
+    tFemtoV0->SetStatusPos(999);
+    tFemtoV0->SetStatusNeg(999);
+  }
   tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
 }
 
@@ -1167,12 +1176,12 @@ AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP
     }
     while (posstack < mcP->GetEntries());
   }
-  
+
   return 0;
 }
 
-void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack, 
-                                                AliFemtoTrack *tFemtoTrack)
+void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
+                                                 AliFemtoTrack *tFemtoTrack)
 {
 
        // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
@@ -1204,25 +1213,30 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
   aodpid[2] = -100000.0;
   aodpid[3] = -100000.0;
   aodpid[4] = -100000.0;
-               
+
   double tTOF = 0.0;
+  Float_t probMis = 1.0;
 
   //what is that code? for what do we need it? nsigma values are not enough?
-   if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) {  //AliESDtrack::kTOFpid=0x8000
-     tTOF = tAodTrack->GetTOFsignal();
-     tAodTrack->GetIntegratedTimes(aodpid);
+  if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) {  //AliESDtrack::kTOFpid=0x8000
+    tTOF = tAodTrack->GetTOFsignal();
+    tAodTrack->GetIntegratedTimes(aodpid);
 
-     tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
-   }
+    tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
+
+    probMis = fAODpidUtil->GetTOFMismatchProbability(tAodTrack);
+  }
+
+
+
+  tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
 
-   tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
   //////  TPC ////////////////////////////////////////////
 
-  float nsigmaTPCK=-1000.;                                                  
-  float nsigmaTPCPi=-1000.;                                                 
-  float nsigmaTPCP=-1000.;                                                  
-          
+  float nsigmaTPCK=-1000.;
+  float nsigmaTPCPi=-1000.;
+  float nsigmaTPCP=-1000.;
+
   //   cout<<"in reader fESDpid"<<fESDpid<<endl;
 
   if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
@@ -1235,26 +1249,27 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
   tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
   tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
 
-  tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());       
-  tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());       
-  tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());      
-  
-  tFemtoTrack->SetTPCsignalN(1); 
-  tFemtoTrack->SetTPCsignalS(1); 
+  tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
+  tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
+  tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
+
+  tFemtoTrack->SetTPCsignalN(1);
+  tFemtoTrack->SetTPCsignalS(1);
   tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
-  ///////TOF//////////////////////
 
-    float vp=-1000.;
-    float nsigmaTOFPi=-1000.;
-    float nsigmaTOFK=-1000.;
-    float nsigmaTOFP=-1000.;
+  ///////TOF//////////////////////
 
-    if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
-       (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
-        (tAodTrack->GetStatus() & AliESDtrack::kTIME) ) //AliESDtrack::kTIME=0x80000000
-      {
-       if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
+  float vp=-1000.;
+  float nsigmaTOFPi=-1000.;
+  float nsigmaTOFK=-1000.;
+  float nsigmaTOFP=-1000.;
+
+  if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
+      (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
+      (tAodTrack->GetStatus() & AliESDtrack::kTIME) //AliESDtrack::kTIME=0x80000000
+      && probMis < 0.01) // TOF mismatch probaility
+  {
+    if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
          {
 
            nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
@@ -1265,21 +1280,21 @@ void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
            Double_t tof=tAodTrack->GetTOFsignal();
            if(tof > 0.) vp=len/tof/0.03;
          }
-      }
-    tFemtoTrack->SetVTOF(vp);
-    tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
-    tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
-    tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
+  }
+  tFemtoTrack->SetVTOF(vp);
+  tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
+  tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
+  tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
+
 
-    
-    //////////////////////////////////////
+  //////////////////////////////////////
 
 }
 
 void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
 {
   fCentRange[0] = min; fCentRange[1] = max;
-  fUsePreCent = 1; 
+  fUsePreCent = 1;
   fEstEventMult = kCentrality;
 }
 
@@ -1313,7 +1328,7 @@ void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s)
 
 void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz)
 {
-    fisEPVZ = iepvz;
+  fisEPVZ = iepvz;
 }
 
 void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
@@ -1364,11 +1379,11 @@ void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrac
 
       // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
       while (globalRadius>Rwanted[iR]){
-       x-=.1;
-       //      printf("propagating to x %5.2f\n",x);
-       if(!etp.PropagateTo(x,bfield))break;
-       etp.GetXYZ(xyz); // GetXYZ returns global coordinates
-       globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
+        x-=.1;
+        //      printf("propagating to x %5.2f\n",x);
+        if(!etp.PropagateTo(x,bfield))break;
+        etp.GetXYZ(xyz); // GetXYZ returns global coordinates
+        globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
       }
       //printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]);
       globalPositionsAtRadii[iR][0]=xyz[0];
index 2f3aa8075070de355ab889e2c4cd67a4a4cdb6b8..f2b64445c4f0fe65e934c0d3c803468cfc6f112b 100644 (file)
@@ -1,13 +1,13 @@
 /*
 ***************************************************************************
 *
-* $Id$ 
+* $Id$
+*
 *
-* 
 ***************************************************************************
 *
-* 
-*              
+*
+*
 *
 ***************************************************************************
 *
@@ -47,7 +47,7 @@
 #include "AliFemtoESDTrackCut.h"
 #include <cstdio>
 
-#ifdef __ROOT__ 
+#ifdef __ROOT__
 ClassImp(AliFemtoESDTrackCut)
 #endif
 
@@ -89,6 +89,8 @@ ClassImp(AliFemtoESDTrackCut)
     fStatus(0),
     fPIDMethod(knSigma),
   fNsigmaTPCTOF(kFALSE),
+  fNsigmaTPConly(kFALSE),
+  fNsigma(3.),
     fminTPCclsF(0),
     fminTPCncls(0),
     fminITScls(0),
@@ -99,7 +101,7 @@ ClassImp(AliFemtoESDTrackCut)
     fNTracksFailed(0),
     fRemoveKinks(kFALSE),
     fRemoveITSFake(kFALSE),
-    fMostProbable(0), 
+    fMostProbable(0),
     fMaxImpactXY(1000.0),
   fMinImpactXY(-1000.0),
     fMaxImpactZ(1000.0),
@@ -132,6 +134,8 @@ ClassImp(AliFemtoESDTrackCut)
   fminITScls=0;
   fPIDMethod=knSigma;
   fNsigmaTPCTOF=kFALSE;
+  fNsigmaTPConly=kFALSE;
+  fNsigma=3.;
 }
 //------------------------------
 AliFemtoESDTrackCut::~AliFemtoESDTrackCut(){
@@ -142,11 +146,11 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
 {
   //cout<<"AliFemtoESDTrackCut::Pass"<<endl;
 
-  // test the particle and return 
+  // test the particle and return
   // true if it meets all the criteria
   // false if it doesn't meet at least one of the criteria
   float tMost[5];
-  
+
   //cout<<"AliFemtoESD  cut"<<endl;
   //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
   if (fStatus!=0)
@@ -157,7 +161,7 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
          //      cout<<track->Flags()<<" "<<fStatus<<" no go through status"<<endl;
          return false;
        }
-       
+
     }
   if (fRemoveKinks) {
     if ((track->KinkIndex(0)) || (track->KinkIndex(1)) || (track->KinkIndex(2)))
@@ -191,12 +195,12 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
 
   if (fMaxImpactZ < TMath::Abs(track->ImpactZ()))
     return false;
-  
+
   if (fMaxSigmaToVertex < track->SigmaToVertex()) {
     return false;
   }
-  
-  if (track->ITSncls() > 0) 
+
+  if (track->ITSncls() > 0)
     if ((track->ITSchi2()/track->ITSncls()) > fMaxITSchiNdof) {
       return false;
     }
@@ -219,13 +223,13 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
          //   cout<<"No Go Through the cut"<<endl;
          //  cout<<fLabel<<" Label="<<track->Label()<<endl;
          return false;
-       }    
+       }
     }
   if (fCharge!=0)
-    {              
+    {
       //cout<<"AliFemtoESD  cut ch "<<endl;
       //cout<<fCharge<<" Charge="<<track->Charge()<<endl;
-      if (track->Charge()!= fCharge)   
+      if (track->Charge()!= fCharge)
        {
          fNTracksFailed++;
          //  cout<<"No Go Through the cut"<<endl;
@@ -235,39 +239,39 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
     }
 
 
-  
+
 
   Bool_t tTPCPidIn = (track->Flags()&AliFemtoTrack::kTPCpid)>0;
   Bool_t tITSPidIn = (track->Flags()&AliFemtoTrack::kITSpid)>0;
   Bool_t tTOFPidIn = (track->Flags()&AliFemtoTrack::kTOFpid)>0;
-  
+
   if(fMinPforTOFpid > 0 && track->P().Mag() > fMinPforTOFpid &&
      track->P().Mag() < fMaxPforTOFpid && !tTOFPidIn)
     {
       fNTracksFailed++;
       return false;
     }
-  
+
   if(fMinPforTPCpid > 0 && track->P().Mag() > fMinPforTPCpid &&
      track->P().Mag() < fMaxPforTPCpid && !tTPCPidIn)
     {
       fNTracksFailed++;
       return false;
     }
-  
+
   if(fMinPforITSpid > 0 && track->P().Mag() > fMinPforITSpid &&
      track->P().Mag() < fMaxPforITSpid && !tITSPidIn)
     {
       fNTracksFailed++;
       return false;
     }
-  
+
 
   float tEnergy = ::sqrt(track->P().Mag2()+fMass*fMass);
   float tRapidity = 0.5*::log((tEnergy+track->P().z())/(tEnergy-track->P().z()));
   float tPt = ::sqrt((track->P().x())*(track->P().x())+(track->P().y())*(track->P().y()));
   float tEta = track->P().PseudoRapidity();
-  
+
   if (fMaxImpactXYPtOff < 999.0) {
     if ((fMaxImpactXYPtOff + fMaxImpactXYPtNrm*TMath::Power(tPt, fMaxImpactXYPtPow)) < TMath::Abs(track->ImpactD())) {
       fNTracksFailed++;
@@ -278,14 +282,14 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
   if ((tRapidity<fRapidity[0])||(tRapidity>fRapidity[1]))
     {
       fNTracksFailed++;
-      //cout<<"No Go Through the cut"<<endl;   
+      //cout<<"No Go Through the cut"<<endl;
       //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
       return false;
     }
   if ((tEta<fEta[0])||(tEta>fEta[1]))
     {
       fNTracksFailed++;
-      //cout<<"No Go Through the cut"<<endl;   
+      //cout<<"No Go Through the cut"<<endl;
       //cout<<fEta[0]<<" < Eta ="<<tEta<<" <"<<fEta[1]<<endl;
       return false;
     }
@@ -300,15 +304,15 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
 
 
 
-  //   cout << "Track has pids: " 
-  //        << track->PidProbElectron() << " " 
-  //        << track->PidProbMuon() << " " 
-  //        << track->PidProbPion() << " " 
-  //        << track->PidProbKaon() << " " 
-  //        << track->PidProbProton() << " " 
+  //   cout << "Track has pids: "
+  //        << track->PidProbElectron() << " "
+  //        << track->PidProbMuon() << " "
+  //        << track->PidProbPion() << " "
+  //        << track->PidProbKaon() << " "
+  //        << track->PidProbProton() << " "
   //        << track->PidProbElectron()+track->PidProbMuon()+track->PidProbPion()+track->PidProbKaon()+track->PidProbProton() << endl;
 
-    
+
   if ((track->PidProbElectron()<fPidProbElectron[0])||(track->PidProbElectron()>fPidProbElectron[1]))
     {
       fNTracksFailed++;
@@ -346,7 +350,7 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
     }
 
   if (fMostProbable) {
-  
+
     int imost=0;
     tMost[0] = track->PidProbElectron()*PidFractionElectron(track->P().Mag());
     tMost[1] = 0.0;
@@ -364,14 +368,14 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
              imost = 2;
 
          }
-         else if (fMostProbable == 3) { 
+         else if (fMostProbable == 3) {
+
 
-    
           if (IsKaonNSigma(track->P().Mag(), track->NSigmaTPCK(), track->NSigmaTOFK())){
-                   
+
              imost = 3;
            }
-  
+
          }
       else if (fMostProbable == 4) { // proton nsigma-PID required contour adjusting (in LHC10h)
         if ( IsProtonNSigma(track->P().Mag(), track->NSigmaTPCP(), track->NSigmaTOFP()) && (TMath::Abs(track->NSigmaTPCP()) < TMath::Abs(track->NSigmaTPCPi())) && (TMath::Abs(track->NSigmaTPCP()) < TMath::Abs(track->NSigmaTPCK())) && (TMath::Abs(track->NSigmaTOFP()) < TMath::Abs(track->NSigmaTOFPi())) && (TMath::Abs(track->NSigmaTOFP()) < TMath::Abs(track->NSigmaTOFK()))
@@ -381,8 +385,8 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
          }
 
        }
-       
-       
+
+
 
     //****Contour Method****
        if(fPIDMethod==1){
@@ -441,14 +445,14 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
                else {
                  if (!(IsKaonTPCdEdx(track->P().Mag(), track->TPCsignal())))
                    imost = 0;
-                 else 
+                 else
                    imost = 3;
                }
              }
            }
            //       }
          }
-    
+
          // Looking for protons
          else if (fMostProbable == 4) {
            //       if (imost == 3) {
@@ -475,7 +479,7 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
                else {
                  if (!(IsKaonTPCdEdx(track->P().Mag(), track->TPCsignal())))
                    imost = 0;
-                 else 
+                 else
                    imost = 3;
                }
              }
@@ -485,7 +489,7 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
        }
     if (imost != fMostProbable) return false;
   }
-  
+
   //fan
   //cout<<"****** Go Through the cut ******"<<endl;
   // cout<<fLabel<<" Label="<<track->Label()<<endl;
@@ -499,8 +503,8 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
   //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
   fNTracksPassed++ ;
   return true;
-    
-    
+
+
 }
 //------------------------------
 AliFemtoString AliFemtoESDTrackCut::Report()
@@ -515,7 +519,7 @@ AliFemtoString AliFemtoESDTrackCut::Report()
   snprintf(tCtemp , 100, "Particle pT:\t%E - %E\n",fPt[0],fPt[1]);
   tStemp+=tCtemp;
   snprintf(tCtemp , 100, "Particle rapidity:\t%E - %E\n",fRapidity[0],fRapidity[1]);
-  tStemp+=tCtemp; 
+  tStemp+=tCtemp;
   snprintf(tCtemp , 100, "Particle eta:\t%E - %E\n",fEta[0],fEta[1]);
   tStemp+=tCtemp;
   snprintf(tCtemp , 100, "Number of tracks which passed:\t%ld  Number which failed:\t%ld\n",fNTracksPassed,fNTracksFailed);
@@ -592,12 +596,12 @@ void AliFemtoESDTrackCut::SetRemoveKinks(const bool& flag)
 {
   fRemoveKinks = flag;
 }
-                           
+
 void AliFemtoESDTrackCut::SetRemoveITSFake(const bool& flag)
 {
   fRemoveITSFake = flag;
 }
-                           
+
 // electron
 // 0.13 - 1.8
 // 0       7.594129e-02    8.256141e-03
@@ -610,30 +614,30 @@ void AliFemtoESDTrackCut::SetRemoveITSFake(const bool& flag)
 float AliFemtoESDTrackCut::PidFractionElectron(float mom) const
 {
   // Provide a parameterized fraction of electrons dependent on momentum
-  if (mom<0.13) 
-    return (7.594129e-02 
-           -5.535827e-01*0.13     
-           +1.728591e+00*0.13*0.13    
-           -2.827893e+00*0.13*0.13*0.13 
-           +2.503553e+00*0.13*0.13*0.13*0.13      
-           -1.125965e+00*0.13*0.13*0.13*0.13*0.13      
-           +2.009036e-01*0.13*0.13*0.13*0.13*0.13*0.13);   
+  if (mom<0.13)
+    return (7.594129e-02
+           -5.535827e-01*0.13
+           +1.728591e+00*0.13*0.13
+           -2.827893e+00*0.13*0.13*0.13
+           +2.503553e+00*0.13*0.13*0.13*0.13
+           -1.125965e+00*0.13*0.13*0.13*0.13*0.13
+           +2.009036e-01*0.13*0.13*0.13*0.13*0.13*0.13);
 
   if (mom>1.8)
-    return (7.594129e-02 
-           -5.535827e-01*1.8      
-           +1.728591e+00*1.8*1.8    
-           -2.827893e+00*1.8*1.8*1.8 
-           +2.503553e+00*1.8*1.8*1.8*1.8          
-           -1.125965e+00*1.8*1.8*1.8*1.8*1.8      
-           +2.009036e-01*1.8*1.8*1.8*1.8*1.8*1.8);   
-  return (7.594129e-02 
-         -5.535827e-01*mom        
-         +1.728591e+00*mom*mom    
-         -2.827893e+00*mom*mom*mom 
-         +2.503553e+00*mom*mom*mom*mom    
-         -1.125965e+00*mom*mom*mom*mom*mom      
-         +2.009036e-01*mom*mom*mom*mom*mom*mom);   
+    return (7.594129e-02
+           -5.535827e-01*1.8
+           +1.728591e+00*1.8*1.8
+           -2.827893e+00*1.8*1.8*1.8
+           +2.503553e+00*1.8*1.8*1.8*1.8
+           -1.125965e+00*1.8*1.8*1.8*1.8*1.8
+           +2.009036e-01*1.8*1.8*1.8*1.8*1.8*1.8);
+  return (7.594129e-02
+         -5.535827e-01*mom
+         +1.728591e+00*mom*mom
+         -2.827893e+00*mom*mom*mom
+         +2.503553e+00*mom*mom*mom*mom
+         -1.125965e+00*mom*mom*mom*mom*mom
+         +2.009036e-01*mom*mom*mom*mom*mom*mom);
 }
 
 // pion
@@ -644,11 +648,11 @@ float AliFemtoESDTrackCut::PidFractionElectron(float mom) const
 float AliFemtoESDTrackCut::PidFractionPion(float mom) const
 {
   // Provide a parameterized fraction of pions dependent on momentum
-  if (mom<0.13) 
+  if (mom<0.13)
     return ( 1.063457e+00
             -4.222208e-01*0.13
             +1.042004e-01*0.0169);
-  if (mom>2.0) 
+  if (mom>2.0)
     return ( 1.063457e+00
             -4.222208e-01*2.0
             +1.042004e-01*4.0);
@@ -666,19 +670,19 @@ float AliFemtoESDTrackCut::PidFractionPion(float mom) const
 float AliFemtoESDTrackCut::PidFractionKaon(float mom) const
 {
   // Provide a parameterized fraction of kaons dependent on momentum
-  if (mom<0.18) 
+  if (mom<0.18)
     return (-7.289406e-02
-           +4.415666e-01*0.18     
-           -2.996790e-01*0.18*0.18    
+           +4.415666e-01*0.18
+           -2.996790e-01*0.18*0.18
            +6.704652e-02*0.18*0.18*0.18);
-  if (mom>2.0) 
+  if (mom>2.0)
     return (-7.289406e-02
-           +4.415666e-01*2.0      
-           -2.996790e-01*2.0*2.0    
+           +4.415666e-01*2.0
+           -2.996790e-01*2.0*2.0
            +6.704652e-02*2.0*2.0*2.0);
   return (-7.289406e-02
-         +4.415666e-01*mom        
-         -2.996790e-01*mom*mom    
+         +4.415666e-01*mom
+         -2.996790e-01*mom*mom
          +6.704652e-02*mom*mom*mom);
 }
 
@@ -692,15 +696,15 @@ float AliFemtoESDTrackCut::PidFractionProton(float mom) const
 {
   // Provide a parameterized fraction of protons dependent on momentum
   if (mom<0.26) return  0.0;
-  if (mom>2.0) 
-    return (-3.730200e-02  
-           +1.163684e-01*2.0         
-           +8.354116e-02*2.0*2.0       
+  if (mom>2.0)
+    return (-3.730200e-02
+           +1.163684e-01*2.0
+           +8.354116e-02*2.0*2.0
            -4.608098e-02*2.0*2.0*2.0);
-  return (-3.730200e-02  
-         +1.163684e-01*mom           
-         +8.354116e-02*mom*mom       
-         -4.608098e-02*mom*mom*mom);  
+  return (-3.730200e-02
+         +1.163684e-01*mom
+         +8.354116e-02*mom*mom
+         -4.608098e-02*mom*mom*mom);
 }
 
 void AliFemtoESDTrackCut::SetMomRangeTOFpidIs(const float& minp, const float& maxp)
@@ -759,12 +763,12 @@ bool AliFemtoESDTrackCut::IsKaonTPCdEdx(float mom, float dEdx)
 //     if (dEdx > a4*mom+b4) return false;
 //     if (dEdx <        b5) return false;
 //   }
-//   else 
+//   else
 //     return false;
 //   //   else {
 //   //     if (dEdx > b5) return false;
 //   //   }
-   
+
 //   return true;
 
   double a1 = -268.896; double b1 =  198.669;
@@ -802,7 +806,7 @@ bool AliFemtoESDTrackCut::IsProtonTPCdEdx(float mom, float dEdx)
   }
 
   return true;
-   
+
 }
 
 bool AliFemtoESDTrackCut::IsPionTOFTime(float mom, float ttof)
@@ -869,7 +873,7 @@ bool AliFemtoESDTrackCut::IsKaonTPCdEdxNSigma(float mom, float nsigmaK)
 
 
   if(mom<0.35 && TMath::Abs(nsigmaK)<5.0)return true;
-  if(mom>=0.35 && mom<0.5 && TMath::Abs(nsigmaK)<3.0)return true; 
+  if(mom>=0.35 && mom<0.5 && TMath::Abs(nsigmaK)<3.0)return true;
   if(mom>=0.5 && mom<0.7 && TMath::Abs(nsigmaK)<2.0)return true;
 
   return false;
@@ -880,7 +884,7 @@ bool AliFemtoESDTrackCut::IsKaonTOFNSigma(float mom, float nsigmaK)
 //  cout<<" AliFemtoESDTrackCut::IsKaonTPCdEdxNSigma "<<mom<<" "<<nsigmaK<<endl;
   //fan
   //  if(mom<1.5 && TMath::Abs(nsigmaK)<3.0)return true;
-  if(mom>=1.5 && TMath::Abs(nsigmaK)<2.0)return true; 
+  if(mom>=1.5 && TMath::Abs(nsigmaK)<2.0)return true;
   return false;
 }
 
@@ -892,19 +896,19 @@ bool AliFemtoESDTrackCut::IsKaonNSigma(float mom, float nsigmaTPCK, float nsigma
   if(mom<0.5)
     {
          if(TMath::Abs(nsigmaTPCK)<2.0)
-          { 
+          {
           return true;
-          } 
-          else 
+          }
+          else
           {
           return false;
           }
     }
-    
-    
+
+
    if(mom>=0.5)
     {
-         if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) 
+         if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0)
          {
          return true;
          }
@@ -913,9 +917,9 @@ bool AliFemtoESDTrackCut::IsKaonNSigma(float mom, float nsigmaTPCK, float nsigma
          return false;
          }
     }
-    
+
 //   if(mom>1.5 || mom<0.15)return false;
-    
+
 
 }
 
@@ -961,7 +965,7 @@ bool AliFemtoESDTrackCut::IsPionNSigma(float mom, float nsigmaTPCPi, float nsigm
          if(mom<0.35 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
          else if(mom<0.5 && mom>=0.35 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
          else if(mom>=0.5 && TMath::Abs(nsigmaTPCPi)<2.0) return true;
-         else return false;      
+         else return false;
        }
       else if(TMath::Abs(nsigmaTOFPi)<3.0 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
     }
@@ -981,14 +985,19 @@ bool AliFemtoESDTrackCut::IsProtonNSigma(float mom, float nsigmaTPCP, float nsig
 
   if (fNsigmaTPCTOF) {
     if (mom > 0.8) {
-        if (TMath::Hypot( nsigmaTOFP, nsigmaTPCP )/TMath::Sqrt(2) < 3.0)
+//        if (TMath::Hypot( nsigmaTOFP, nsigmaTPCP )/TMath::Sqrt(2) < 3.0)
+        if (TMath::Hypot( nsigmaTOFP, nsigmaTPCP ) < fNsigma)
             return true;
        }
     else {
-        if (TMath::Abs(nsigmaTPCP) < 3.0)
+        if (TMath::Abs(nsigmaTPCP) < fNsigma)
             return true;
     }
   }
+  else if (fNsigmaTPConly) {
+    if (TMath::Abs(nsigmaTPCP) < fNsigma)
+      return true;
+  }
   else {
     if (mom > 0.8 && mom < 2.5) {
       if ( TMath::Abs(nsigmaTPCP) < 3.0 && TMath::Abs(nsigmaTOFP) < 3.0)
@@ -1018,16 +1027,26 @@ void AliFemtoESDTrackCut::SetNsigmaTPCTOF(Bool_t nsigma)
 {
   fNsigmaTPCTOF = nsigma;
 }
-void AliFemtoESDTrackCut::SetClusterRequirementITS(AliESDtrackCuts::Detector det, AliESDtrackCuts::ITSClusterRequirement req) 
-{ 
-  fCutClusterRequirementITS[det] = req; 
+
+void AliFemtoESDTrackCut::SetNsigmaTPConly(Bool_t nsigma)
+{
+  fNsigmaTPConly = nsigma;
+}
+
+void AliFemtoESDTrackCut::SetNsigma(Double_t nsigma)
+{
+  fNsigma = nsigma;
+}
+
+void AliFemtoESDTrackCut::SetClusterRequirementITS(AliESDtrackCuts::Detector det, AliESDtrackCuts::ITSClusterRequirement req)
+{
+  fCutClusterRequirementITS[det] = req;
 }
 
 Bool_t AliFemtoESDTrackCut::CheckITSClusterRequirement(AliESDtrackCuts::ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
 {
   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
-  
+
   switch (req)
     {
     case AliESDtrackCuts::kOff:        return kTRUE;
@@ -1039,6 +1058,6 @@ Bool_t AliFemtoESDTrackCut::CheckITSClusterRequirement(AliESDtrackCuts::ITSClust
     case AliESDtrackCuts::kOnlySecond: return clusterL2 && !clusterL1;
     case AliESDtrackCuts::kBoth:       return clusterL1 && clusterL2;
   }
-  
+
   return kFALSE;
 }
index 5252547d2769ec307955f5515e9306212b5940d4..9518ec2d322bfc09d023792cd063318f1627895a 100644 (file)
@@ -1,7 +1,7 @@
 ///////////////////////////////////////////////////////////////////////////
 //                                                                       //
 // AliFemtoESDTrackCut: A basic track cut that used information from     //
-// ALICE ESD to accept or reject the track.                              //  
+// ALICE ESD to accept or reject the track.                              //
 // Enables the selection on charge, transverse momentum, rapidity,       //
 // pid probabilities, number of ITS and TPC clusters                     //
 // Author: Marek Chojnacki (WUT), mchojnacki@knf.pw.edu.pl               //
 #include "AliFemtoTrackCut.h"
 
 
-class AliFemtoESDTrackCut : public AliFemtoTrackCut 
+class AliFemtoESDTrackCut : public AliFemtoTrackCut
 {
   public:
 
   enum PIDMethodType {knSigma=0, kContour=1};
-  typedef enum PIDMethodType ReadPIDMethodType; 
+  typedef enum PIDMethodType ReadPIDMethodType;
 
   AliFemtoESDTrackCut();
   virtual ~AliFemtoESDTrackCut();
@@ -62,9 +62,11 @@ class AliFemtoESDTrackCut : public AliFemtoTrackCut
   void SetMostProbablePion();
   void SetMostProbableKaon();
   void SetMostProbableProton();
-  void SetNoMostProbable(); 
+  void SetNoMostProbable();
   void SetPIDMethod(ReadPIDMethodType newMethod);
   void SetNsigmaTPCTOF(Bool_t);
+  void SetNsigmaTPConly(Bool_t);
+  void SetNsigma(Double_t);
   void SetClusterRequirementITS(AliESDtrackCuts::Detector det, AliESDtrackCuts::ITSClusterRequirement req = AliESDtrackCuts::kOff);
 
   void SetMomRangeTOFpidIs(const float& minp, const float& maxp);
@@ -81,17 +83,19 @@ class AliFemtoESDTrackCut : public AliFemtoTrackCut
   float             fPidProbPion[2];     // bounds for pion probability
   float             fPidProbKaon[2];     // bounds for kaon probability
   float             fPidProbProton[2];   // bounds for proton probability
-  float             fPidProbMuon[2];     // bounds for muon probability 
+  float             fPidProbMuon[2];     // bounds for muon probability
 
   AliESDtrackCuts::ITSClusterRequirement fCutClusterRequirementITS[3];  // detailed ITS cluster requirements for (SPD, SDD, SSD) - from AliESDtrackcuts!
-  bool              fLabel;              // if true label<0 will not pass throught 
+  bool              fLabel;              // if true label<0 will not pass throught
   long              fStatus;             // staus flag
-  ReadPIDMethodType fPIDMethod;          // which PID mehod to use. 0 - nsgima, 1 - contour 
+  ReadPIDMethodType fPIDMethod;          // which PID mehod to use. 0 - nsgima, 1 - contour
   Bool_t            fNsigmaTPCTOF;       // true if squared nsigma from TPC and TOF, false if separately from TPC and TOF
+  Bool_t            fNsigmaTPConly;       // true if nsigma from TPC only
+  Double_t            fNsigma;       // number of sigmas - 3 by default
 
   short             fminTPCclsF;         // min number of findable clusters in the TPC
   short             fminTPCncls;         // min number of clusters in the TPC
-  int               fminITScls;          // min number of clusters assigned in the ITS 
+  int               fminITScls;          // min number of clusters assigned in the ITS
   float             fMaxITSchiNdof;      // maximum allowed chi2/ndof for ITS clusters
   float             fMaxTPCchiNdof;      // maximum allowed chi2/ndof for TPC clusters
   float             fMaxSigmaToVertex;   // maximum allowed sigma to primary vertex
@@ -138,7 +142,7 @@ class AliFemtoESDTrackCut : public AliFemtoTrackCut
   Bool_t CheckITSClusterRequirement(AliESDtrackCuts::ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2); //the same as in AliESDtrackCuts
 
 
-#ifdef __ROOT__ 
+#ifdef __ROOT__
   ClassDef(AliFemtoESDTrackCut, 1)
 #endif
     };
diff --git a/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/AvgSep/ConfigFemtoAnalysis.C b/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/AvgSep/ConfigFemtoAnalysis.C
new file mode 100644 (file)
index 0000000..b7540f8
--- /dev/null
@@ -0,0 +1,449 @@
+/*********************************************************************
+ *                                                                   *
+ * Configfemtoanalysis.C - configuration macro for the femtoscopic   *
+ * analysis, meant as a QA process for two-particle effects          *
+ *                                                                   *
+ * Author: Adam Kisiel (Adam.Kisiel@cern.ch)                         *
+ *                                                                   *
+ *********************************************************************/
+
+#if !defined(__CINT__) || defined(__MAKECINT_)
+#include "AliFemtoManager.h"
+#include "AliFemtoEventReaderESDChain.h"
+#include "AliFemtoEventReaderESDChainKine.h"
+#include "AliFemtoEventReaderAODChain.h"
+#include "AliFemtoSimpleAnalysis.h"
+#include "AliFemtoBasicEventCut.h"
+#include "AliFemtoESDTrackCut.h"
+#include "AliFemtoCorrFctn.h"
+#include "AliFemtoCutMonitorParticleYPt.h"
+#include "AliFemtoCutMonitorParticleVertPos.h"
+#include "AliFemtoCutMonitorParticleMomRes.h"
+#include "AliFemtoCutMonitorParticlePID.h"
+#include "AliFemtoCutMonitorEventMult.h"
+#include "AliFemtoCutMonitorEventVertex.h"
+#include "AliFemtoShareQualityTPCEntranceSepPairCut.h"
+#include "AliFemtoPairCutAntiGamma.h"
+#include "AliFemtoPairCutRadialDistance.h"
+#include "AliFemtoQinvCorrFctn.h"
+#include "AliFemtoCorrFctnNonIdDR.h"
+#include "AliFemtoShareQualityCorrFctn.h"
+#include "AliFemtoTPCInnerCorrFctn.h"
+#include "AliFemtoVertexMultAnalysis.h"
+#include "AliFemtoCorrFctn3DSpherical.h"
+#include "AliFemtoChi2CorrFctn.h"
+#include "AliFemtoCorrFctnTPCNcls.h"
+#include "AliFemtoBPLCMS3DCorrFctn.h"
+#include "AliFemtoCorrFctn3DLCMSSym.h"
+#include "AliFemtoModelBPLCMSCorrFctn.h"
+#include "AliFemtoModelCorrFctn3DSpherical.h"
+#include "AliFemtoModelGausLCMSFreezeOutGenerator.h"
+#include "AliFemtoModelGausRinvFreezeOutGenerator.h"
+#include "AliFemtoModelManager.h"
+#include "AliFemtoModelWeightGeneratorBasic.h"
+#include "AliFemtoModelWeightGeneratorLednicky.h"
+#include "AliFemtoCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnSource.h"
+#include "AliFemtoCutMonitorParticlePtPDG.h"
+#include "AliFemtoKTPairCut.h"
+#include "AliFemtoAvgSepCorrFctn.h"
+#endif
+
+//________________________________________________________________________
+AliFemtoManager* ConfigFemtoAnalysis() {
+
+  double PionMass = 0.13956995;
+  double KaonMass = 0.493677;
+  double ProtonMass = 0.938272013;
+
+  // double psi = TMath::Pi()/2.;
+  // double psid = TMath::Pi()/6.;
+
+  // int runepvzero[7] = {1, 1, 1, 1, 1, 1, 1};
+  // double epvzerobins[7] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid, -psi+4*psid, -psi+5*psid, -psi+6*psid};
+
+  double psi = TMath::Pi()/2.;
+  double psid = TMath::Pi()/3.;
+
+  int runepvzero[4] = {0, 0, 0, 1};
+  double epvzerobins[4] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid};
+
+  int runmults[10] = {1, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  int multbins[11] = {0.001, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900};
+
+  int runch[3] = {1, 0, 1};
+  const char *chrgs[3] = { "PP", "APAP", "PAP" };
+
+  int runktdep = 0;
+  double ktrng[3] = {0.01, 1.0, 5.0};
+
+  int numOfMultBins = 10;
+  int numOfChTypes = 3;
+  int numOfkTbins = 2;
+  int numOfEPvzero = 4;
+
+  int runqinv = 1;
+  int runshlcms = 0;// 0:PRF(PAP), 1:LCMS(PP,APAP)
+
+  int runtype = 2; // Types 0 - global, 1 - ITS only, 2 - TPC Inner
+  int isrealdata = 1;
+
+  //  int gammacut = 1;
+
+  double shqmax = 1.0;
+  int nbinssh = 100;
+
+  AliFemtoEventReaderAODChain *Reader = new AliFemtoEventReaderAODChain();
+  Reader->SetFilterBit(7);
+  Reader->SetCentralityPreSelection(0.001, 510);
+  Reader->SetEPVZERO(kTRUE);
+
+  AliFemtoManager* Manager = new AliFemtoManager();
+  Manager->SetEventReader(Reader);
+
+  AliFemtoVertexMultAnalysis    *anetaphitpc[10*3*2];
+  AliFemtoBasicEventCut         *mecetaphitpc[10*3*2];
+  AliFemtoCutMonitorEventMult   *cutPassEvMetaphitpc[50];
+  AliFemtoCutMonitorEventMult   *cutFailEvMetaphitpc[50];
+  // AliFemtoCutMonitorEventVertex *cutPassEvVetaphitpc[50];
+  // AliFemtoCutMonitorEventVertex *cutFailEvVetaphitpc[50];
+  AliFemtoESDTrackCut           *dtc1etaphitpc[50];
+  AliFemtoESDTrackCut           *dtc2etaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass2PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail2PIDetaphitpc[50];
+  //   AliFemtoPairCutAntiGamma      *sqpcetaphitpcdiff[10*3];
+  //   AliFemtoShareQualityTPCEntranceSepPairCut      *sqpcetaphitpcsame[10*3];
+  //AliFemtoPairCutAntiGamma      *sqpcetaphitpc[10*3];
+  AliFemtoPairCutRadialDistance      *sqpcetaphitpc[50];
+  //  AliFemtoChi2CorrFctn          *cchiqinvetaphitpc[20*2];
+  AliFemtoKTPairCut             *ktpcuts[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmtpc[50];
+  AliFemtoCorrFctnDirectYlm     *cylmkttpc[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmetaphitpc[10*3];
+  AliFemtoQinvCorrFctn          *cqinvkttpc[50*2];
+  AliFemtoQinvCorrFctn          *cqinvtpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstartpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstarkttpc[50*2];
+  AliFemtoCorrFctnDEtaDPhi      *cdedpetaphi[50*2];
+  AliFemtoAvgSepCorrFctn          *cAvgSeptpc[50];
+
+  //   AliFemtoCorrFctn3DLCMSSym     *cq3dlcmskttpc[20*2];
+  //   AliFemtoCorrFctnTPCNcls       *cqinvnclstpc[20];
+  //   AliFemtoShareQualityCorrFctn  *cqinvsqtpc[20*10];
+  //   AliFemtoChi2CorrFctn          *cqinvchi2tpc[20];
+    AliFemtoTPCInnerCorrFctn      *cqinvinnertpc[50];
+
+  // *** Third QA task - HBT analysis with all pair cuts off, TPC only ***
+  // *** Begin pion-pion (positive) analysis ***
+  int aniter = 0;
+
+  for (int imult = 0; imult < numOfMultBins; imult++) {
+    if (runmults[imult]) {
+
+      for (int ichg = 0; ichg < numOfChTypes; ichg++) {
+        if (runch[ichg]) {
+
+          for (int iepvzero = 0; iepvzero < numOfEPvzero; iepvzero++) {
+            if (runepvzero[iepvzero]) {
+
+              aniter = imult * numOfChTypes + ichg * numOfEPvzero + iepvzero;
+              // aniter = ichg * numOfMultBins + imult * numOfEPvzero + iepvzero;
+
+              // cout << "aniter = " << aniter << endl;
+              //        aniter = ichg * numOfMultBins + imult;
+
+              // if (ichg == 2)
+              //       runshlcms = 0;
+              // else
+              //       runshlcms = 1;
+
+
+              //________________________
+
+              anetaphitpc[aniter] = new AliFemtoVertexMultAnalysis(8, -8.0, 8.0, 4, multbins[imult], multbins[imult+1]);
+              anetaphitpc[aniter]->SetNumEventsToMix(10);
+              anetaphitpc[aniter]->SetMinSizePartCollection(1);
+              anetaphitpc[aniter]->SetVerboseMode(kFALSE);
+
+              mecetaphitpc[aniter] = new AliFemtoBasicEventCut();
+              mecetaphitpc[aniter]->SetEventMult(0.001,100000);
+              mecetaphitpc[aniter]->SetVertZPos(-8,8);
+
+              if (iepvzero == 3)
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[0],epvzerobins[3]);
+              else
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[iepvzero],epvzerobins[iepvzero+1]);
+
+              // if (isrealdata)
+              //     mecetaphitpc[aniter]->SetAcceptOnlyPhysics(kTRUE);
+
+              // cutPassEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutPass%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              // cutFailEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutFail%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              // mecetaphitpc[aniter]->AddCutMonitor(cutPassEvMetaphitpc[aniter], cutFailEvMetaphitpc[aniter]);
+
+              // cutPassEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutPass%stpcM%i", chrgs[ichg], imult));
+              // cutFailEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutFail%stpcM%i", chrgs[ichg], imult));
+              // mecetaphitpc[aniter]->AddCutMonitor(cutPassEvVetaphitpc[aniter], cutFailEvVetaphitpc[aniter]);
+
+              dtc1etaphitpc[aniter] = new AliFemtoESDTrackCut();
+              dtc2etaphitpc[aniter] = new AliFemtoESDTrackCut();
+
+              if (ichg == 0) {
+                dtc1etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 1) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 2) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc2etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+                dtc2etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+
+              dtc1etaphitpc[aniter]->SetEta(-0.8,0.8);
+              dtc1etaphitpc[aniter]->SetMass(ProtonMass);
+              dtc1etaphitpc[aniter]->SetMostProbableProton();
+              dtc1etaphitpc[aniter]->SetNsigma(3.0);
+              //dtc1etaphitpc[aniter]->SetNsigma(2.0);
+              dtc1etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+              //dtc1etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              if (ichg == 2) {
+                dtc2etaphitpc[aniter]->SetEta(-0.8,0.8);
+                dtc2etaphitpc[aniter]->SetMass(ProtonMass);
+                dtc2etaphitpc[aniter]->SetMostProbableProton();
+                dtc2etaphitpc[aniter]->SetNsigma(3.0);
+                //dtc2etaphitpc[aniter]->SetNsigma(2.0);
+                dtc2etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+                //dtc2etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              }
+
+              // Track quality cuts
+
+              if (runtype == 0) {
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                //            dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0182, 0.0350, -1.01);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 1) {
+                //      dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit|AliESDtrack::kITSpureSA);
+                //      dtc1etaphitpc[aniter]->SetminTPCncls(70);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                //      dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.25);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 2) {
+                //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                //dtc1etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                // dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+
+                if (ichg == 2) {
+                  //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                  dtc2etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                  dtc2etaphitpc[aniter]->SetminTPCncls(80);
+                  dtc2etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                  dtc2etaphitpc[aniter]->SetLabel(kFALSE);
+                  dtc2etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                  //dtc2etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                  // dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+                }
+
+              }
+
+              // cutPass1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              // cutFail1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              // dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1YPtetaphitpc[aniter], cutFail1YPtetaphitpc[aniter]);
+
+              // cutPass1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2);//0-pion,1-kaon,2-proton
+              // cutFail1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult , iepvzero),2);
+              // dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1PIDetaphitpc[aniter], cutFail1PIDetaphitpc[aniter]);
+
+              // if (ichg == 2){
+              //     cutPass2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass2%stpcM%i", chrgs[ichg], imult),2);//0-pion,1-kaon,2-proton
+              //     cutFail2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail2%stpcM%i", chrgs[ichg], imult),2);
+              //     dtc2etaphitpc[aniter]->AddCutMonitor(cutPass2PIDetaphitpc[aniter], cutFail2PIDetaphitpc[aniter]);
+              // }
+
+              // sqpcetaphitpc[aniter] = new AliFemtoPairCutAntiGamma();
+              sqpcetaphitpc[aniter] = new AliFemtoPairCutRadialDistance();
+
+              if (runtype == 0) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(1.5);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(0.12, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 1) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(1.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(5.0);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(1.2, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 2) {
+                //sqpcetaphitpc[aniter]->SetUseAOD(kTRUE);
+
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+
+                //         if (gammacut == 0) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //}
+                //else if (gammacut == 1) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //}
+
+                // sqpcetaphitpc[aniter]->SetMagneticFieldSign(-1); // field1 -1, field3 +1
+                // sqpcetaphitpc[aniter]->SetMinimumRadius(0.8); // biggest inefficiency for R=1.1 m (checked on small sample)
+
+                // // sqpcetaphitpc[aniter]->SetMinimumRadius(1.6); //0.8
+                // sqpcetaphitpc[aniter]->SetPhiStarMin(kTRUE);
+                // sqpcetaphitpc[aniter]->SetPhiStarDifferenceMinimum(0.04); // 0.012 - pions, 0.017 - kaons, 0.018
+                // sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.03); // 0.017 - pions, 0.015 - kaons
+
+              }
+
+              anetaphitpc[aniter]->SetEventCut(mecetaphitpc[aniter]);
+
+              if (ichg == 2) {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc2etaphitpc[aniter]);
+              }
+              else {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc1etaphitpc[aniter]);
+              }
+
+              anetaphitpc[aniter]->SetPairCut(sqpcetaphitpc[aniter]);
+
+
+              // if (ichg == 2) {
+              //   ckstartpc[aniter] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),nbinssh,0.0,shqmax);
+              //   anetaphitpc[aniter]->AddCorrFctn(ckstartpc[aniter]);
+              // }
+              // else {
+
+              //   cqinvtpc[aniter] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2*nbinssh,0.0,2*shqmax);
+              //   anetaphitpc[aniter]->AddCorrFctn(cqinvtpc[aniter]);
+
+              // }
+
+              // cylmtpc[aniter] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%i", chrgs[ichg], imult),2,nbinssh, 0.0,shqmax,runshlcms);
+              // anetaphitpc[aniter]->AddCorrFctn(cylmtpc[aniter]);
+
+
+              cAvgSeptpc[aniter] = new AliFemtoAvgSepCorrFctn(Form("cAvgSep%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),4*nbinssh,0.0,200);
+              anetaphitpc[aniter]->AddCorrFctn(cAvgSeptpc[aniter]);
+
+              cqinvinnertpc[aniter] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%d", chrgs[ichg], imult),nbinssh,0.0,shqmax);
+              cqinvinnertpc[aniter]->SetRadius(1.2);
+              anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[aniter]);
+
+
+              if (runktdep) {
+                int ktm;
+                for (int ikt=0; ikt<numOfkTbins; ikt++) {
+
+                  ktm = aniter * numOfkTbins + ikt;
+                  ktpcuts[ktm] = new AliFemtoKTPairCut(ktrng[ikt], ktrng[ikt+1]);
+
+
+                  // cylmkttpc[ktm] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%ikT%i", chrgs[ichg], imult, ikt),2,nbinssh,0.0,shqmax,runshlcms);
+                  // cylmkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  // anetaphitpc[aniter]->AddCorrFctn(cylmkttpc[ktm]);
+
+
+
+                  if (ichg == 2) {
+                    ckstarkttpc[ktm] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),nbinssh,0.0,shqmax);
+                    ckstarkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(ckstarkttpc[ktm]);
+
+
+                  }
+                  else {
+                    cqinvkttpc[ktm] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),2*nbinssh,0.0,2*shqmax);
+                    cqinvkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(cqinvkttpc[ktm]);
+                  }
+
+                  //         cqinvsqtpc[ktm] = new AliFemtoShareQualityCorrFctn(Form("cqinvsq%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  //         cqinvsqtpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //         anetaphitpc[aniter]->AddCorrFctn(cqinvsqtpc[ktm]);
+
+                  // cqinvinnertpc[ktm] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  // cqinvinnertpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  // cqinvinnertpc[ktm]->SetRadius(1.2);
+                  // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[ktm]);
+
+                  //         if (run3d) {
+                  //           cq3dlcmskttpc[ktm] = new AliFemtoCorrFctn3DLCMSSym(Form("cq3d%stpcM%ikT%i", chrgs[ichg], imult, ikt),60,(imult>3)?((imult>6)?((imult>7)?0.6:0.4):0.25):0.15);
+                  //           cq3dlcmskttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //           anetaphitpc[aniter]->AddCorrFctn(cq3dlcmskttpc[ktm]);
+                  //         }
+                }
+              }
+
+              // cdedpetaphi[aniter] = new AliFemtoCorrFctnDEtaDPhi(Form("cdedp%stpcM%i", chrgs[ichg], imult),240, 240);
+              // anetaphitpc[aniter]->AddCorrFctn(cdedpetaphi[aniter]);
+
+              Manager->AddAnalysis(anetaphitpc[aniter]);
+            }
+          }
+        }
+      }
+    }
+  }
+
+
+// *** End pion-pion (positive) analysis
+
+
+  return Manager;
+}
diff --git a/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcOnly/ConfigFemtoAnalysis.C b/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcOnly/ConfigFemtoAnalysis.C
new file mode 100644 (file)
index 0000000..a7ea8fa
--- /dev/null
@@ -0,0 +1,447 @@
+/*********************************************************************
+ *                                                                   *
+ * Configfemtoanalysis.C - configuration macro for the femtoscopic   *
+ * analysis, meant as a QA process for two-particle effects          *
+ *                                                                   *
+ * Author: Adam Kisiel (Adam.Kisiel@cern.ch)                         *
+ *                                                                   *
+ *********************************************************************/
+
+#if !defined(__CINT__) || defined(__MAKECINT_)
+#include "AliFemtoManager.h"
+#include "AliFemtoEventReaderESDChain.h"
+#include "AliFemtoEventReaderESDChainKine.h"
+#include "AliFemtoEventReaderAODChain.h"
+#include "AliFemtoSimpleAnalysis.h"
+#include "AliFemtoBasicEventCut.h"
+#include "AliFemtoESDTrackCut.h"
+#include "AliFemtoCorrFctn.h"
+#include "AliFemtoCutMonitorParticleYPt.h"
+#include "AliFemtoCutMonitorParticleVertPos.h"
+#include "AliFemtoCutMonitorParticleMomRes.h"
+#include "AliFemtoCutMonitorParticlePID.h"
+#include "AliFemtoCutMonitorEventMult.h"
+#include "AliFemtoCutMonitorEventVertex.h"
+#include "AliFemtoShareQualityTPCEntranceSepPairCut.h"
+#include "AliFemtoPairCutAntiGamma.h"
+#include "AliFemtoPairCutRadialDistance.h"
+#include "AliFemtoQinvCorrFctn.h"
+#include "AliFemtoCorrFctnNonIdDR.h"
+#include "AliFemtoShareQualityCorrFctn.h"
+#include "AliFemtoTPCInnerCorrFctn.h"
+#include "AliFemtoVertexMultAnalysis.h"
+#include "AliFemtoCorrFctn3DSpherical.h"
+#include "AliFemtoChi2CorrFctn.h"
+#include "AliFemtoCorrFctnTPCNcls.h"
+#include "AliFemtoBPLCMS3DCorrFctn.h"
+#include "AliFemtoCorrFctn3DLCMSSym.h"
+#include "AliFemtoModelBPLCMSCorrFctn.h"
+#include "AliFemtoModelCorrFctn3DSpherical.h"
+#include "AliFemtoModelGausLCMSFreezeOutGenerator.h"
+#include "AliFemtoModelGausRinvFreezeOutGenerator.h"
+#include "AliFemtoModelManager.h"
+#include "AliFemtoModelWeightGeneratorBasic.h"
+#include "AliFemtoModelWeightGeneratorLednicky.h"
+#include "AliFemtoCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnSource.h"
+#include "AliFemtoCutMonitorParticlePtPDG.h"
+#include "AliFemtoKTPairCut.h"
+#include "AliFemtoAvgSepCorrFctn.h"
+#endif
+
+//________________________________________________________________________
+AliFemtoManager* ConfigFemtoAnalysis() {
+
+  double PionMass = 0.13956995;
+  double KaonMass = 0.493677;
+  double ProtonMass = 0.938272013;
+
+  // double psi = TMath::Pi()/2.;
+  // double psid = TMath::Pi()/6.;
+
+  // int runepvzero[7] = {1, 1, 1, 1, 1, 1, 1};
+  // double epvzerobins[7] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid, -psi+4*psid, -psi+5*psid, -psi+6*psid};
+
+  double psi = TMath::Pi()/2.;
+  double psid = TMath::Pi()/3.;
+
+  int runepvzero[4] = {0, 0, 0, 1};
+  double epvzerobins[4] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid};
+
+  int runmults[10] = {1, 1, 0, 0, 0, 0, 0, 0, 0, 0};
+  int multbins[11] = {0.001, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900};
+
+  int runch[3] = {1, 1, 1};
+  const char *chrgs[3] = { "PP", "APAP", "PAP" };
+
+  int runktdep = 1;
+  double ktrng[3] = {0.01, 1.0, 5.0};
+
+  int numOfMultBins = 10;
+  int numOfChTypes = 3;
+  int numOfkTbins = 2;
+  int numOfEPvzero = 4;
+
+  int runqinv = 1;
+  int runshlcms = 0;// 0:PRF(PAP), 1:LCMS(PP,APAP)
+
+  int runtype = 2; // Types 0 - global, 1 - ITS only, 2 - TPC Inner
+  int isrealdata = 1;
+
+  //  int gammacut = 1;
+
+  double shqmax = 1.0;
+  int nbinssh = 100;
+
+  AliFemtoEventReaderAODChain *Reader = new AliFemtoEventReaderAODChain();
+  Reader->SetFilterBit(7);
+  Reader->SetCentralityPreSelection(0.001, 510);
+  Reader->SetEPVZERO(kTRUE);
+
+  AliFemtoManager* Manager = new AliFemtoManager();
+  Manager->SetEventReader(Reader);
+
+  AliFemtoVertexMultAnalysis    *anetaphitpc[10*3*2];
+  AliFemtoBasicEventCut         *mecetaphitpc[10*3*2];
+  AliFemtoCutMonitorEventMult   *cutPassEvMetaphitpc[50];
+  AliFemtoCutMonitorEventMult   *cutFailEvMetaphitpc[50];
+  // AliFemtoCutMonitorEventVertex *cutPassEvVetaphitpc[50];
+  // AliFemtoCutMonitorEventVertex *cutFailEvVetaphitpc[50];
+  AliFemtoESDTrackCut           *dtc1etaphitpc[50];
+  AliFemtoESDTrackCut           *dtc2etaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass2PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail2PIDetaphitpc[50];
+  //   AliFemtoPairCutAntiGamma      *sqpcetaphitpcdiff[10*3];
+  //   AliFemtoShareQualityTPCEntranceSepPairCut      *sqpcetaphitpcsame[10*3];
+  //AliFemtoPairCutAntiGamma      *sqpcetaphitpc[10*3];
+  AliFemtoPairCutRadialDistance      *sqpcetaphitpc[50];
+  //  AliFemtoChi2CorrFctn          *cchiqinvetaphitpc[20*2];
+  AliFemtoKTPairCut             *ktpcuts[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmtpc[50];
+  AliFemtoCorrFctnDirectYlm     *cylmkttpc[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmetaphitpc[10*3];
+  AliFemtoQinvCorrFctn          *cqinvkttpc[50*2];
+  AliFemtoQinvCorrFctn          *cqinvtpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstartpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstarkttpc[50*2];
+  AliFemtoCorrFctnDEtaDPhi      *cdedpetaphi[50*2];
+  AliFemtoAvgSepCorrFctn          *cAvgSeptpc[50];
+
+  //   AliFemtoCorrFctn3DLCMSSym     *cq3dlcmskttpc[20*2];
+  //   AliFemtoCorrFctnTPCNcls       *cqinvnclstpc[20];
+  //   AliFemtoShareQualityCorrFctn  *cqinvsqtpc[20*10];
+  //   AliFemtoChi2CorrFctn          *cqinvchi2tpc[20];
+    AliFemtoTPCInnerCorrFctn      *cqinvinnertpc[50];
+
+  // *** Third QA task - HBT analysis with all pair cuts off, TPC only ***
+  // *** Begin pion-pion (positive) analysis ***
+  int aniter = 0;
+
+  for (int imult = 0; imult < numOfMultBins; imult++) {
+    if (runmults[imult]) {
+
+      for (int ichg = 0; ichg < numOfChTypes; ichg++) {
+        if (runch[ichg]) {
+
+          for (int iepvzero = 0; iepvzero < numOfEPvzero; iepvzero++) {
+            if (runepvzero[iepvzero]) {
+
+              aniter = imult * numOfChTypes + ichg * numOfEPvzero + iepvzero;
+              // aniter = ichg * numOfMultBins + imult * numOfEPvzero + iepvzero;
+
+              // cout << "aniter = " << aniter << endl;
+              //        aniter = ichg * numOfMultBins + imult;
+
+              // if (ichg == 2)
+              //       runshlcms = 0;
+              // else
+              //       runshlcms = 1;
+
+
+              //________________________
+
+              anetaphitpc[aniter] = new AliFemtoVertexMultAnalysis(8, -8.0, 8.0, 4, multbins[imult], multbins[imult+1]);
+              anetaphitpc[aniter]->SetNumEventsToMix(10);
+              anetaphitpc[aniter]->SetMinSizePartCollection(1);
+              anetaphitpc[aniter]->SetVerboseMode(kFALSE);
+
+              mecetaphitpc[aniter] = new AliFemtoBasicEventCut();
+              mecetaphitpc[aniter]->SetEventMult(0.001,100000);
+              mecetaphitpc[aniter]->SetVertZPos(-8,8);
+
+              if (iepvzero == 3)
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[0],epvzerobins[3]);
+              else
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[iepvzero],epvzerobins[iepvzero+1]);
+
+              // if (isrealdata)
+              //     mecetaphitpc[aniter]->SetAcceptOnlyPhysics(kTRUE);
+
+              // cutPassEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutPass%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              // cutFailEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutFail%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              // mecetaphitpc[aniter]->AddCutMonitor(cutPassEvMetaphitpc[aniter], cutFailEvMetaphitpc[aniter]);
+
+              // cutPassEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutPass%stpcM%i", chrgs[ichg], imult));
+              // cutFailEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutFail%stpcM%i", chrgs[ichg], imult));
+              // mecetaphitpc[aniter]->AddCutMonitor(cutPassEvVetaphitpc[aniter], cutFailEvVetaphitpc[aniter]);
+
+              dtc1etaphitpc[aniter] = new AliFemtoESDTrackCut();
+              dtc2etaphitpc[aniter] = new AliFemtoESDTrackCut();
+
+              if (ichg == 0) {
+                dtc1etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 1) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 2) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc2etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+                dtc2etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+
+              dtc1etaphitpc[aniter]->SetEta(-0.8,0.8);
+              dtc1etaphitpc[aniter]->SetMass(ProtonMass);
+              dtc1etaphitpc[aniter]->SetMostProbableProton();
+              dtc1etaphitpc[aniter]->SetNsigma(3.0);
+              //dtc1etaphitpc[aniter]->SetNsigma(2.0);
+              //dtc1etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+              dtc1etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              if (ichg == 2) {
+                dtc2etaphitpc[aniter]->SetEta(-0.8,0.8);
+                dtc2etaphitpc[aniter]->SetMass(ProtonMass);
+                dtc2etaphitpc[aniter]->SetMostProbableProton();
+                dtc2etaphitpc[aniter]->SetNsigma(3.0);
+                //dtc2etaphitpc[aniter]->SetNsigma(2.0);
+                //dtc2etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+                dtc2etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              }
+
+              // Track quality cuts
+
+              if (runtype == 0) {
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                //            dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0182, 0.0350, -1.01);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 1) {
+                //      dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit|AliESDtrack::kITSpureSA);
+                //      dtc1etaphitpc[aniter]->SetminTPCncls(70);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                //      dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.25);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 2) {
+                //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                //dtc1etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                // dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+
+                if (ichg == 2) {
+                  //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                  dtc2etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                  dtc2etaphitpc[aniter]->SetminTPCncls(80);
+                  dtc2etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                  dtc2etaphitpc[aniter]->SetLabel(kFALSE);
+                  dtc2etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                  //dtc2etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                  // dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+                }
+
+              }
+
+              // cutPass1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              // cutFail1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              // dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1YPtetaphitpc[aniter], cutFail1YPtetaphitpc[aniter]);
+
+              // cutPass1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2);//0-pion,1-kaon,2-proton
+              // cutFail1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult , iepvzero),2);
+              // dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1PIDetaphitpc[aniter], cutFail1PIDetaphitpc[aniter]);
+
+              // if (ichg == 2){
+              //     cutPass2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass2%stpcM%i", chrgs[ichg], imult),2);//0-pion,1-kaon,2-proton
+              //     cutFail2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail2%stpcM%i", chrgs[ichg], imult),2);
+              //     dtc2etaphitpc[aniter]->AddCutMonitor(cutPass2PIDetaphitpc[aniter], cutFail2PIDetaphitpc[aniter]);
+              // }
+
+              // sqpcetaphitpc[aniter] = new AliFemtoPairCutAntiGamma();
+              sqpcetaphitpc[aniter] = new AliFemtoPairCutRadialDistance();
+
+              if (runtype == 0) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(1.5);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(0.12, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 1) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(1.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(5.0);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(1.2, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 2) {
+                //sqpcetaphitpc[aniter]->SetUseAOD(kTRUE);
+
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+
+                //         if (gammacut == 0) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //}
+                //else if (gammacut == 1) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //}
+
+                // sqpcetaphitpc[aniter]->SetMagneticFieldSign(-1); // field1 -1, field3 +1
+                // sqpcetaphitpc[aniter]->SetMinimumRadius(0.8); // biggest inefficiency for R=1.1 m (checked on small sample)
+
+                sqpcetaphitpc[aniter]->SetMinimumRadius(1.2); //0.8
+                sqpcetaphitpc[aniter]->SetPhiStarMin(kFALSE);
+                sqpcetaphitpc[aniter]->SetPhiStarDifferenceMinimum(0.045); // 0.012 - pions, 0.017 - kaons, 0.018
+                sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.01); // 0.017 - pions, 0.015 - kaons
+
+              }
+
+              anetaphitpc[aniter]->SetEventCut(mecetaphitpc[aniter]);
+
+              if (ichg == 2) {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc2etaphitpc[aniter]);
+              }
+              else {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc1etaphitpc[aniter]);
+              }
+
+              anetaphitpc[aniter]->SetPairCut(sqpcetaphitpc[aniter]);
+
+
+              if (ichg == 2) {
+                ckstartpc[aniter] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),nbinssh,0.0,shqmax);
+                anetaphitpc[aniter]->AddCorrFctn(ckstartpc[aniter]);
+              }
+              else {
+
+                cqinvtpc[aniter] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2*nbinssh,0.0,2*shqmax);
+                anetaphitpc[aniter]->AddCorrFctn(cqinvtpc[aniter]);
+
+              }
+
+              cylmtpc[aniter] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%i", chrgs[ichg], imult),2,nbinssh, 0.0,shqmax,runshlcms);
+              anetaphitpc[aniter]->AddCorrFctn(cylmtpc[aniter]);
+
+
+              // cAvgSeptpc[aniter] = new AliFemtoAvgSepCorrFctn(Form("cAvgSep%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),4*nbinssh,0.0,200);
+              // anetaphitpc[aniter]->AddCorrFctn(cAvgSeptpc[aniter]);
+
+              // cqinvinnertpc[aniter] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%d", chrgs[ichg], imult),nbinssh,0.0,shqmax);
+              // cqinvinnertpc[aniter]->SetRadius(1.2);
+              // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[aniter]);
+
+
+              if (runktdep) {
+                int ktm;
+                for (int ikt=0; ikt<numOfkTbins; ikt++) {
+
+                  ktm = aniter * numOfkTbins + ikt;
+                  ktpcuts[ktm] = new AliFemtoKTPairCut(ktrng[ikt], ktrng[ikt+1]);
+
+
+                  cylmkttpc[ktm] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%ikT%i", chrgs[ichg], imult, ikt),2,nbinssh,0.0,shqmax,runshlcms);
+                  cylmkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  anetaphitpc[aniter]->AddCorrFctn(cylmkttpc[ktm]);
+
+                  if (ichg == 2) {
+                    ckstarkttpc[ktm] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),nbinssh,0.0,shqmax);
+                    ckstarkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(ckstarkttpc[ktm]);
+
+
+                  }
+                  else {
+                    cqinvkttpc[ktm] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),2*nbinssh,0.0,2*shqmax);
+                    cqinvkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(cqinvkttpc[ktm]);
+                  }
+
+                  //         cqinvsqtpc[ktm] = new AliFemtoShareQualityCorrFctn(Form("cqinvsq%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  //         cqinvsqtpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //         anetaphitpc[aniter]->AddCorrFctn(cqinvsqtpc[ktm]);
+
+                  // cqinvinnertpc[ktm] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  // cqinvinnertpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  // cqinvinnertpc[ktm]->SetRadius(1.2);
+                  // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[ktm]);
+
+                  //         if (run3d) {
+                  //           cq3dlcmskttpc[ktm] = new AliFemtoCorrFctn3DLCMSSym(Form("cq3d%stpcM%ikT%i", chrgs[ichg], imult, ikt),60,(imult>3)?((imult>6)?((imult>7)?0.6:0.4):0.25):0.15);
+                  //           cq3dlcmskttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //           anetaphitpc[aniter]->AddCorrFctn(cq3dlcmskttpc[ktm]);
+                  //         }
+                }
+              }
+
+              // cdedpetaphi[aniter] = new AliFemtoCorrFctnDEtaDPhi(Form("cdedp%stpcM%i", chrgs[ichg], imult),240, 240);
+              // anetaphitpc[aniter]->AddCorrFctn(cdedpetaphi[aniter]);
+
+              Manager->AddAnalysis(anetaphitpc[aniter]);
+            }
+          }
+        }
+      }
+    }
+  }
+
+
+// *** End pion-pion (positive) analysis
+
+
+  return Manager;
+}
diff --git a/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcOnly/monitors/ConfigFemtoAnalysis.C b/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcOnly/monitors/ConfigFemtoAnalysis.C
new file mode 100644 (file)
index 0000000..d3df0a9
--- /dev/null
@@ -0,0 +1,447 @@
+/*********************************************************************
+ *                                                                   *
+ * Configfemtoanalysis.C - configuration macro for the femtoscopic   *
+ * analysis, meant as a QA process for two-particle effects          *
+ *                                                                   *
+ * Author: Adam Kisiel (Adam.Kisiel@cern.ch)                         *
+ *                                                                   *
+ *********************************************************************/
+
+#if !defined(__CINT__) || defined(__MAKECINT_)
+#include "AliFemtoManager.h"
+#include "AliFemtoEventReaderESDChain.h"
+#include "AliFemtoEventReaderESDChainKine.h"
+#include "AliFemtoEventReaderAODChain.h"
+#include "AliFemtoSimpleAnalysis.h"
+#include "AliFemtoBasicEventCut.h"
+#include "AliFemtoESDTrackCut.h"
+#include "AliFemtoCorrFctn.h"
+#include "AliFemtoCutMonitorParticleYPt.h"
+#include "AliFemtoCutMonitorParticleVertPos.h"
+#include "AliFemtoCutMonitorParticleMomRes.h"
+#include "AliFemtoCutMonitorParticlePID.h"
+#include "AliFemtoCutMonitorEventMult.h"
+#include "AliFemtoCutMonitorEventVertex.h"
+#include "AliFemtoShareQualityTPCEntranceSepPairCut.h"
+#include "AliFemtoPairCutAntiGamma.h"
+#include "AliFemtoPairCutRadialDistance.h"
+#include "AliFemtoQinvCorrFctn.h"
+#include "AliFemtoCorrFctnNonIdDR.h"
+#include "AliFemtoShareQualityCorrFctn.h"
+#include "AliFemtoTPCInnerCorrFctn.h"
+#include "AliFemtoVertexMultAnalysis.h"
+#include "AliFemtoCorrFctn3DSpherical.h"
+#include "AliFemtoChi2CorrFctn.h"
+#include "AliFemtoCorrFctnTPCNcls.h"
+#include "AliFemtoBPLCMS3DCorrFctn.h"
+#include "AliFemtoCorrFctn3DLCMSSym.h"
+#include "AliFemtoModelBPLCMSCorrFctn.h"
+#include "AliFemtoModelCorrFctn3DSpherical.h"
+#include "AliFemtoModelGausLCMSFreezeOutGenerator.h"
+#include "AliFemtoModelGausRinvFreezeOutGenerator.h"
+#include "AliFemtoModelManager.h"
+#include "AliFemtoModelWeightGeneratorBasic.h"
+#include "AliFemtoModelWeightGeneratorLednicky.h"
+#include "AliFemtoCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnSource.h"
+#include "AliFemtoCutMonitorParticlePtPDG.h"
+#include "AliFemtoKTPairCut.h"
+#include "AliFemtoAvgSepCorrFctn.h"
+#endif
+
+//________________________________________________________________________
+AliFemtoManager* ConfigFemtoAnalysis() {
+
+  double PionMass = 0.13956995;
+  double KaonMass = 0.493677;
+  double ProtonMass = 0.938272013;
+
+  // double psi = TMath::Pi()/2.;
+  // double psid = TMath::Pi()/6.;
+
+  // int runepvzero[7] = {1, 1, 1, 1, 1, 1, 1};
+  // double epvzerobins[7] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid, -psi+4*psid, -psi+5*psid, -psi+6*psid};
+
+  double psi = TMath::Pi()/2.;
+  double psid = TMath::Pi()/3.;
+
+  int runepvzero[4] = {0, 0, 0, 1};
+  double epvzerobins[4] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid};
+
+  int runmults[10] = {1, 1, 0, 0, 0, 0, 0, 0, 0, 0};
+  int multbins[11] = {0.001, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900};
+
+  int runch[3] = {1, 1, 1};
+  const char *chrgs[3] = { "PP", "APAP", "PAP" };
+
+  int runktdep = 0;
+  double ktrng[3] = {0.01, 1.0, 5.0};
+
+  int numOfMultBins = 10;
+  int numOfChTypes = 3;
+  int numOfkTbins = 2;
+  int numOfEPvzero = 4;
+
+  int runqinv = 1;
+  int runshlcms = 0;// 0:PRF(PAP), 1:LCMS(PP,APAP)
+
+  int runtype = 2; // Types 0 - global, 1 - ITS only, 2 - TPC Inner
+  int isrealdata = 1;
+
+  //  int gammacut = 1;
+
+  double shqmax = 1.0;
+  int nbinssh = 100;
+
+  AliFemtoEventReaderAODChain *Reader = new AliFemtoEventReaderAODChain();
+  Reader->SetFilterBit(7);
+  Reader->SetCentralityPreSelection(0.001, 510);
+  Reader->SetEPVZERO(kTRUE);
+
+  AliFemtoManager* Manager = new AliFemtoManager();
+  Manager->SetEventReader(Reader);
+
+  AliFemtoVertexMultAnalysis    *anetaphitpc[10*3*2];
+  AliFemtoBasicEventCut         *mecetaphitpc[10*3*2];
+  AliFemtoCutMonitorEventMult   *cutPassEvMetaphitpc[50];
+  AliFemtoCutMonitorEventMult   *cutFailEvMetaphitpc[50];
+  AliFemtoCutMonitorEventVertex *cutPassEvVetaphitpc[50];
+   AliFemtoCutMonitorEventVertex *cutFailEvVetaphitpc[50];
+  AliFemtoESDTrackCut           *dtc1etaphitpc[50];
+  AliFemtoESDTrackCut           *dtc2etaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass2PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail2PIDetaphitpc[50];
+  //   AliFemtoPairCutAntiGamma      *sqpcetaphitpcdiff[10*3];
+  //   AliFemtoShareQualityTPCEntranceSepPairCut      *sqpcetaphitpcsame[10*3];
+  //AliFemtoPairCutAntiGamma      *sqpcetaphitpc[10*3];
+  AliFemtoPairCutRadialDistance      *sqpcetaphitpc[50];
+  //  AliFemtoChi2CorrFctn          *cchiqinvetaphitpc[20*2];
+  AliFemtoKTPairCut             *ktpcuts[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmtpc[50];
+  AliFemtoCorrFctnDirectYlm     *cylmkttpc[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmetaphitpc[10*3];
+  AliFemtoQinvCorrFctn          *cqinvkttpc[50*2];
+  AliFemtoQinvCorrFctn          *cqinvtpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstartpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstarkttpc[50*2];
+  AliFemtoCorrFctnDEtaDPhi      *cdedpetaphi[50*2];
+  AliFemtoAvgSepCorrFctn          *cAvgSeptpc[50];
+
+  //   AliFemtoCorrFctn3DLCMSSym     *cq3dlcmskttpc[20*2];
+  //   AliFemtoCorrFctnTPCNcls       *cqinvnclstpc[20];
+  //   AliFemtoShareQualityCorrFctn  *cqinvsqtpc[20*10];
+  //   AliFemtoChi2CorrFctn          *cqinvchi2tpc[20];
+    AliFemtoTPCInnerCorrFctn      *cqinvinnertpc[50];
+
+  // *** Third QA task - HBT analysis with all pair cuts off, TPC only ***
+  // *** Begin pion-pion (positive) analysis ***
+  int aniter = 0;
+
+  for (int imult = 0; imult < numOfMultBins; imult++) {
+    if (runmults[imult]) {
+
+      for (int ichg = 0; ichg < numOfChTypes; ichg++) {
+        if (runch[ichg]) {
+
+          for (int iepvzero = 0; iepvzero < numOfEPvzero; iepvzero++) {
+            if (runepvzero[iepvzero]) {
+
+              aniter = imult * numOfChTypes + ichg * numOfEPvzero + iepvzero;
+              // aniter = ichg * numOfMultBins + imult * numOfEPvzero + iepvzero;
+
+              // cout << "aniter = " << aniter << endl;
+              //        aniter = ichg * numOfMultBins + imult;
+
+              // if (ichg == 2)
+              //       runshlcms = 0;
+              // else
+              //       runshlcms = 1;
+
+
+              //________________________
+
+              anetaphitpc[aniter] = new AliFemtoVertexMultAnalysis(8, -8.0, 8.0, 4, multbins[imult], multbins[imult+1]);
+              anetaphitpc[aniter]->SetNumEventsToMix(10);
+              anetaphitpc[aniter]->SetMinSizePartCollection(1);
+              anetaphitpc[aniter]->SetVerboseMode(kFALSE);
+
+              mecetaphitpc[aniter] = new AliFemtoBasicEventCut();
+              mecetaphitpc[aniter]->SetEventMult(0.001,100000);
+              mecetaphitpc[aniter]->SetVertZPos(-8,8);
+
+              if (iepvzero == 3)
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[0],epvzerobins[3]);
+              else
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[iepvzero],epvzerobins[iepvzero+1]);
+
+              // if (isrealdata)
+              //     mecetaphitpc[aniter]->SetAcceptOnlyPhysics(kTRUE);
+
+              cutPassEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutPass%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              cutFailEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutFail%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              mecetaphitpc[aniter]->AddCutMonitor(cutPassEvMetaphitpc[aniter], cutFailEvMetaphitpc[aniter]);
+
+              cutPassEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutPass%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              cutFailEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutFail%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              mecetaphitpc[aniter]->AddCutMonitor(cutPassEvVetaphitpc[aniter], cutFailEvVetaphitpc[aniter]);
+
+              dtc1etaphitpc[aniter] = new AliFemtoESDTrackCut();
+              dtc2etaphitpc[aniter] = new AliFemtoESDTrackCut();
+
+              if (ichg == 0) {
+                dtc1etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 1) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 2) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc2etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+                dtc2etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+
+              dtc1etaphitpc[aniter]->SetEta(-0.8,0.8);
+              dtc1etaphitpc[aniter]->SetMass(ProtonMass);
+              dtc1etaphitpc[aniter]->SetMostProbableProton();
+              dtc1etaphitpc[aniter]->SetNsigma(3.0);
+              //dtc1etaphitpc[aniter]->SetNsigma(2.0);
+              //dtc1etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+              dtc1etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              if (ichg == 2) {
+                dtc2etaphitpc[aniter]->SetEta(-0.8,0.8);
+                dtc2etaphitpc[aniter]->SetMass(ProtonMass);
+                dtc2etaphitpc[aniter]->SetMostProbableProton();
+                dtc2etaphitpc[aniter]->SetNsigma(3.0);
+                //dtc2etaphitpc[aniter]->SetNsigma(2.0);
+                //dtc2etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+                dtc2etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              }
+
+              // Track quality cuts
+
+              if (runtype == 0) {
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                //            dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0182, 0.0350, -1.01);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 1) {
+                //      dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit|AliESDtrack::kITSpureSA);
+                //      dtc1etaphitpc[aniter]->SetminTPCncls(70);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                //      dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.25);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 2) {
+                //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                //dtc1etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                // dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+
+                if (ichg == 2) {
+                  //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                  dtc2etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                  dtc2etaphitpc[aniter]->SetminTPCncls(80);
+                  dtc2etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                  dtc2etaphitpc[aniter]->SetLabel(kFALSE);
+                  dtc2etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                  //dtc2etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                  // dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+                }
+
+              }
+
+              cutPass1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              cutFail1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1YPtetaphitpc[aniter], cutFail1YPtetaphitpc[aniter]);
+
+              cutPass1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2);//0-pion,1-kaon,2-proton
+              cutFail1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult , iepvzero),2);
+              dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1PIDetaphitpc[aniter], cutFail1PIDetaphitpc[aniter]);
+
+              // if (ichg == 2){
+              //     cutPass2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass2%stpcM%i", chrgs[ichg], imult),2);//0-pion,1-kaon,2-proton
+              //     cutFail2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail2%stpcM%i", chrgs[ichg], imult),2);
+              //     dtc2etaphitpc[aniter]->AddCutMonitor(cutPass2PIDetaphitpc[aniter], cutFail2PIDetaphitpc[aniter]);
+              // }
+
+              // sqpcetaphitpc[aniter] = new AliFemtoPairCutAntiGamma();
+              sqpcetaphitpc[aniter] = new AliFemtoPairCutRadialDistance();
+
+              if (runtype == 0) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(1.5);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(0.12, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 1) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(1.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(5.0);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(1.2, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 2) {
+                //sqpcetaphitpc[aniter]->SetUseAOD(kTRUE);
+
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+
+                //         if (gammacut == 0) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //}
+                //else if (gammacut == 1) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //}
+
+                // sqpcetaphitpc[aniter]->SetMagneticFieldSign(-1); // field1 -1, field3 +1
+                // sqpcetaphitpc[aniter]->SetMinimumRadius(0.8); // biggest inefficiency for R=1.1 m (checked on small sample)
+
+                sqpcetaphitpc[aniter]->SetMinimumRadius(1.2); //0.8
+                sqpcetaphitpc[aniter]->SetPhiStarMin(kFALSE);
+                sqpcetaphitpc[aniter]->SetPhiStarDifferenceMinimum(0.045); // 0.012 - pions, 0.017 - kaons, 0.018
+                sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.01); // 0.017 - pions, 0.015 - kaons
+
+              }
+
+              anetaphitpc[aniter]->SetEventCut(mecetaphitpc[aniter]);
+
+              if (ichg == 2) {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc2etaphitpc[aniter]);
+              }
+              else {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc1etaphitpc[aniter]);
+              }
+
+              anetaphitpc[aniter]->SetPairCut(sqpcetaphitpc[aniter]);
+
+
+              // if (ichg == 2) {
+              //   ckstartpc[aniter] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),nbinssh,0.0,shqmax);
+              //   anetaphitpc[aniter]->AddCorrFctn(ckstartpc[aniter]);
+              // }
+              // else {
+
+              //   cqinvtpc[aniter] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2*nbinssh,0.0,2*shqmax);
+              //   anetaphitpc[aniter]->AddCorrFctn(cqinvtpc[aniter]);
+
+              // }
+
+              // cylmtpc[aniter] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%i", chrgs[ichg], imult),2,nbinssh, 0.0,shqmax,runshlcms);
+              // anetaphitpc[aniter]->AddCorrFctn(cylmtpc[aniter]);
+
+
+              // cAvgSeptpc[aniter] = new AliFemtoAvgSepCorrFctn(Form("cAvgSep%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),4*nbinssh,0.0,200);
+              // anetaphitpc[aniter]->AddCorrFctn(cAvgSeptpc[aniter]);
+
+              // cqinvinnertpc[aniter] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%d", chrgs[ichg], imult),nbinssh,0.0,shqmax);
+              // cqinvinnertpc[aniter]->SetRadius(1.2);
+              // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[aniter]);
+
+
+              if (runktdep) {
+                int ktm;
+                for (int ikt=0; ikt<numOfkTbins; ikt++) {
+
+                  ktm = aniter * numOfkTbins + ikt;
+                  ktpcuts[ktm] = new AliFemtoKTPairCut(ktrng[ikt], ktrng[ikt+1]);
+
+
+                  cylmkttpc[ktm] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%ikT%i", chrgs[ichg], imult, ikt),2,nbinssh,0.0,shqmax,runshlcms);
+                  cylmkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  anetaphitpc[aniter]->AddCorrFctn(cylmkttpc[ktm]);
+
+                  if (ichg == 2) {
+                    ckstarkttpc[ktm] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),nbinssh,0.0,shqmax);
+                    ckstarkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(ckstarkttpc[ktm]);
+
+
+                  }
+                  else {
+                    cqinvkttpc[ktm] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),2*nbinssh,0.0,2*shqmax);
+                    cqinvkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(cqinvkttpc[ktm]);
+                  }
+
+                  //         cqinvsqtpc[ktm] = new AliFemtoShareQualityCorrFctn(Form("cqinvsq%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  //         cqinvsqtpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //         anetaphitpc[aniter]->AddCorrFctn(cqinvsqtpc[ktm]);
+
+                  // cqinvinnertpc[ktm] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  // cqinvinnertpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  // cqinvinnertpc[ktm]->SetRadius(1.2);
+                  // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[ktm]);
+
+                  //         if (run3d) {
+                  //           cq3dlcmskttpc[ktm] = new AliFemtoCorrFctn3DLCMSSym(Form("cq3d%stpcM%ikT%i", chrgs[ichg], imult, ikt),60,(imult>3)?((imult>6)?((imult>7)?0.6:0.4):0.25):0.15);
+                  //           cq3dlcmskttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //           anetaphitpc[aniter]->AddCorrFctn(cq3dlcmskttpc[ktm]);
+                  //         }
+                }
+              }
+
+              // cdedpetaphi[aniter] = new AliFemtoCorrFctnDEtaDPhi(Form("cdedp%stpcM%i", chrgs[ichg], imult),240, 240);
+              // anetaphitpc[aniter]->AddCorrFctn(cdedpetaphi[aniter]);
+
+              Manager->AddAnalysis(anetaphitpc[aniter]);
+            }
+          }
+        }
+      }
+    }
+  }
+
+
+// *** End pion-pion (positive) analysis
+
+
+  return Manager;
+}
diff --git a/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcTof/central/ConfigFemtoAnalysis.C b/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcTof/central/ConfigFemtoAnalysis.C
new file mode 100644 (file)
index 0000000..6c8baeb
--- /dev/null
@@ -0,0 +1,447 @@
+/*********************************************************************
+ *                                                                   *
+ * Configfemtoanalysis.C - configuration macro for the femtoscopic   *
+ * analysis, meant as a QA process for two-particle effects          *
+ *                                                                   *
+ * Author: Adam Kisiel (Adam.Kisiel@cern.ch)                         *
+ *                                                                   *
+ *********************************************************************/
+
+#if !defined(__CINT__) || defined(__MAKECINT_)
+#include "AliFemtoManager.h"
+#include "AliFemtoEventReaderESDChain.h"
+#include "AliFemtoEventReaderESDChainKine.h"
+#include "AliFemtoEventReaderAODChain.h"
+#include "AliFemtoSimpleAnalysis.h"
+#include "AliFemtoBasicEventCut.h"
+#include "AliFemtoESDTrackCut.h"
+#include "AliFemtoCorrFctn.h"
+#include "AliFemtoCutMonitorParticleYPt.h"
+#include "AliFemtoCutMonitorParticleVertPos.h"
+#include "AliFemtoCutMonitorParticleMomRes.h"
+#include "AliFemtoCutMonitorParticlePID.h"
+#include "AliFemtoCutMonitorEventMult.h"
+#include "AliFemtoCutMonitorEventVertex.h"
+#include "AliFemtoShareQualityTPCEntranceSepPairCut.h"
+#include "AliFemtoPairCutAntiGamma.h"
+#include "AliFemtoPairCutRadialDistance.h"
+#include "AliFemtoQinvCorrFctn.h"
+#include "AliFemtoCorrFctnNonIdDR.h"
+#include "AliFemtoShareQualityCorrFctn.h"
+#include "AliFemtoTPCInnerCorrFctn.h"
+#include "AliFemtoVertexMultAnalysis.h"
+#include "AliFemtoCorrFctn3DSpherical.h"
+#include "AliFemtoChi2CorrFctn.h"
+#include "AliFemtoCorrFctnTPCNcls.h"
+#include "AliFemtoBPLCMS3DCorrFctn.h"
+#include "AliFemtoCorrFctn3DLCMSSym.h"
+#include "AliFemtoModelBPLCMSCorrFctn.h"
+#include "AliFemtoModelCorrFctn3DSpherical.h"
+#include "AliFemtoModelGausLCMSFreezeOutGenerator.h"
+#include "AliFemtoModelGausRinvFreezeOutGenerator.h"
+#include "AliFemtoModelManager.h"
+#include "AliFemtoModelWeightGeneratorBasic.h"
+#include "AliFemtoModelWeightGeneratorLednicky.h"
+#include "AliFemtoCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnSource.h"
+#include "AliFemtoCutMonitorParticlePtPDG.h"
+#include "AliFemtoKTPairCut.h"
+#include "AliFemtoAvgSepCorrFctn.h"
+#endif
+
+//________________________________________________________________________
+AliFemtoManager* ConfigFemtoAnalysis() {
+
+  double PionMass = 0.13956995;
+  double KaonMass = 0.493677;
+  double ProtonMass = 0.938272013;
+
+  // double psi = TMath::Pi()/2.;
+  // double psid = TMath::Pi()/6.;
+
+  // int runepvzero[7] = {1, 1, 1, 1, 1, 1, 1};
+  // double epvzerobins[7] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid, -psi+4*psid, -psi+5*psid, -psi+6*psid};
+
+  double psi = TMath::Pi()/2.;
+  double psid = TMath::Pi()/3.;
+
+  int runepvzero[4] = {0, 0, 0, 1};
+  double epvzerobins[4] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid};
+
+  int runmults[10] = {1, 1, 0, 0, 0, 0, 0, 0, 0, 0};
+  int multbins[11] = {0.001, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900};
+
+  int runch[3] = {1, 1, 1};
+  const char *chrgs[3] = { "PP", "APAP", "PAP" };
+
+  int runktdep = 1;
+  double ktrng[3] = {0.01, 1.0, 5.0};
+
+  int numOfMultBins = 10;
+  int numOfChTypes = 3;
+  int numOfkTbins = 2;
+  int numOfEPvzero = 4;
+
+  int runqinv = 1;
+  int runshlcms = 0;// 0:PRF(PAP), 1:LCMS(PP,APAP)
+
+  int runtype = 2; // Types 0 - global, 1 - ITS only, 2 - TPC Inner
+  int isrealdata = 1;
+
+  //  int gammacut = 1;
+
+  double shqmax = 1.0;
+  int nbinssh = 100;
+
+  AliFemtoEventReaderAODChain *Reader = new AliFemtoEventReaderAODChain();
+  Reader->SetFilterBit(7);
+  Reader->SetCentralityPreSelection(0.001, 510);
+  Reader->SetEPVZERO(kTRUE);
+
+  AliFemtoManager* Manager = new AliFemtoManager();
+  Manager->SetEventReader(Reader);
+
+  AliFemtoVertexMultAnalysis    *anetaphitpc[10*3*2];
+  AliFemtoBasicEventCut         *mecetaphitpc[10*3*2];
+  AliFemtoCutMonitorEventMult   *cutPassEvMetaphitpc[50];
+  AliFemtoCutMonitorEventMult   *cutFailEvMetaphitpc[50];
+  // AliFemtoCutMonitorEventVertex *cutPassEvVetaphitpc[50];
+  // AliFemtoCutMonitorEventVertex *cutFailEvVetaphitpc[50];
+  AliFemtoESDTrackCut           *dtc1etaphitpc[50];
+  AliFemtoESDTrackCut           *dtc2etaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass2PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail2PIDetaphitpc[50];
+  //   AliFemtoPairCutAntiGamma      *sqpcetaphitpcdiff[10*3];
+  //   AliFemtoShareQualityTPCEntranceSepPairCut      *sqpcetaphitpcsame[10*3];
+  //AliFemtoPairCutAntiGamma      *sqpcetaphitpc[10*3];
+  AliFemtoPairCutRadialDistance      *sqpcetaphitpc[50];
+  //  AliFemtoChi2CorrFctn          *cchiqinvetaphitpc[20*2];
+  AliFemtoKTPairCut             *ktpcuts[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmtpc[50];
+  AliFemtoCorrFctnDirectYlm     *cylmkttpc[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmetaphitpc[10*3];
+  AliFemtoQinvCorrFctn          *cqinvkttpc[50*2];
+  AliFemtoQinvCorrFctn          *cqinvtpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstartpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstarkttpc[50*2];
+  AliFemtoCorrFctnDEtaDPhi      *cdedpetaphi[50*2];
+  AliFemtoAvgSepCorrFctn          *cAvgSeptpc[50];
+
+  //   AliFemtoCorrFctn3DLCMSSym     *cq3dlcmskttpc[20*2];
+  //   AliFemtoCorrFctnTPCNcls       *cqinvnclstpc[20];
+  //   AliFemtoShareQualityCorrFctn  *cqinvsqtpc[20*10];
+  //   AliFemtoChi2CorrFctn          *cqinvchi2tpc[20];
+    AliFemtoTPCInnerCorrFctn      *cqinvinnertpc[50];
+
+  // *** Third QA task - HBT analysis with all pair cuts off, TPC only ***
+  // *** Begin pion-pion (positive) analysis ***
+  int aniter = 0;
+
+  for (int imult = 0; imult < numOfMultBins; imult++) {
+    if (runmults[imult]) {
+
+      for (int ichg = 0; ichg < numOfChTypes; ichg++) {
+        if (runch[ichg]) {
+
+          for (int iepvzero = 0; iepvzero < numOfEPvzero; iepvzero++) {
+            if (runepvzero[iepvzero]) {
+
+              aniter = imult * numOfChTypes + ichg * numOfEPvzero + iepvzero;
+              // aniter = ichg * numOfMultBins + imult * numOfEPvzero + iepvzero;
+
+              // cout << "aniter = " << aniter << endl;
+              //        aniter = ichg * numOfMultBins + imult;
+
+              // if (ichg == 2)
+              //       runshlcms = 0;
+              // else
+              //       runshlcms = 1;
+
+
+              //________________________
+
+              anetaphitpc[aniter] = new AliFemtoVertexMultAnalysis(8, -8.0, 8.0, 4, multbins[imult], multbins[imult+1]);
+              anetaphitpc[aniter]->SetNumEventsToMix(10);
+              anetaphitpc[aniter]->SetMinSizePartCollection(1);
+              anetaphitpc[aniter]->SetVerboseMode(kFALSE);
+
+              mecetaphitpc[aniter] = new AliFemtoBasicEventCut();
+              mecetaphitpc[aniter]->SetEventMult(0.001,100000);
+              mecetaphitpc[aniter]->SetVertZPos(-8,8);
+
+              if (iepvzero == 3)
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[0],epvzerobins[3]);
+              else
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[iepvzero],epvzerobins[iepvzero+1]);
+
+              // if (isrealdata)
+              //     mecetaphitpc[aniter]->SetAcceptOnlyPhysics(kTRUE);
+
+              // cutPassEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutPass%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              // cutFailEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutFail%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              // mecetaphitpc[aniter]->AddCutMonitor(cutPassEvMetaphitpc[aniter], cutFailEvMetaphitpc[aniter]);
+
+              // cutPassEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutPass%stpcM%i", chrgs[ichg], imult));
+              // cutFailEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutFail%stpcM%i", chrgs[ichg], imult));
+              // mecetaphitpc[aniter]->AddCutMonitor(cutPassEvVetaphitpc[aniter], cutFailEvVetaphitpc[aniter]);
+
+              dtc1etaphitpc[aniter] = new AliFemtoESDTrackCut();
+              dtc2etaphitpc[aniter] = new AliFemtoESDTrackCut();
+
+              if (ichg == 0) {
+                dtc1etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 1) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 2) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc2etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+                dtc2etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+
+              dtc1etaphitpc[aniter]->SetEta(-0.8,0.8);
+              dtc1etaphitpc[aniter]->SetMass(ProtonMass);
+              dtc1etaphitpc[aniter]->SetMostProbableProton();
+              dtc1etaphitpc[aniter]->SetNsigma(3.0);
+              //dtc1etaphitpc[aniter]->SetNsigma(2.0);
+              dtc1etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+              //dtc1etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              if (ichg == 2) {
+                dtc2etaphitpc[aniter]->SetEta(-0.8,0.8);
+                dtc2etaphitpc[aniter]->SetMass(ProtonMass);
+                dtc2etaphitpc[aniter]->SetMostProbableProton();
+                dtc2etaphitpc[aniter]->SetNsigma(3.0);
+                //dtc2etaphitpc[aniter]->SetNsigma(2.0);
+                dtc2etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+                //dtc2etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              }
+
+              // Track quality cuts
+
+              if (runtype == 0) {
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                //            dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0182, 0.0350, -1.01);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 1) {
+                //      dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit|AliESDtrack::kITSpureSA);
+                //      dtc1etaphitpc[aniter]->SetminTPCncls(70);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                //      dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.25);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 2) {
+                //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                //dtc1etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                // dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+
+                if (ichg == 2) {
+                  //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                  dtc2etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                  dtc2etaphitpc[aniter]->SetminTPCncls(80);
+                  dtc2etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                  dtc2etaphitpc[aniter]->SetLabel(kFALSE);
+                  dtc2etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                  //dtc2etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                  // dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+                }
+
+              }
+
+              // cutPass1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              // cutFail1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              // dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1YPtetaphitpc[aniter], cutFail1YPtetaphitpc[aniter]);
+
+              // cutPass1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2);//0-pion,1-kaon,2-proton
+              // cutFail1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult , iepvzero),2);
+              // dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1PIDetaphitpc[aniter], cutFail1PIDetaphitpc[aniter]);
+
+              // if (ichg == 2){
+              //     cutPass2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass2%stpcM%i", chrgs[ichg], imult),2);//0-pion,1-kaon,2-proton
+              //     cutFail2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail2%stpcM%i", chrgs[ichg], imult),2);
+              //     dtc2etaphitpc[aniter]->AddCutMonitor(cutPass2PIDetaphitpc[aniter], cutFail2PIDetaphitpc[aniter]);
+              // }
+
+              // sqpcetaphitpc[aniter] = new AliFemtoPairCutAntiGamma();
+              sqpcetaphitpc[aniter] = new AliFemtoPairCutRadialDistance();
+
+              if (runtype == 0) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(1.5);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(0.12, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 1) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(1.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(5.0);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(1.2, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 2) {
+                //sqpcetaphitpc[aniter]->SetUseAOD(kTRUE);
+
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+
+                //         if (gammacut == 0) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //}
+                //else if (gammacut == 1) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //}
+
+                // sqpcetaphitpc[aniter]->SetMagneticFieldSign(-1); // field1 -1, field3 +1
+                // sqpcetaphitpc[aniter]->SetMinimumRadius(0.8); // biggest inefficiency for R=1.1 m (checked on small sample)
+
+                sqpcetaphitpc[aniter]->SetMinimumRadius(1.2); //0.8
+                sqpcetaphitpc[aniter]->SetPhiStarMin(kFALSE);
+                sqpcetaphitpc[aniter]->SetPhiStarDifferenceMinimum(0.045); // 0.012 - pions, 0.017 - kaons, 0.018
+                sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.01); // 0.017 - pions, 0.015 - kaons
+
+              }
+
+              anetaphitpc[aniter]->SetEventCut(mecetaphitpc[aniter]);
+
+              if (ichg == 2) {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc2etaphitpc[aniter]);
+              }
+              else {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc1etaphitpc[aniter]);
+              }
+
+              anetaphitpc[aniter]->SetPairCut(sqpcetaphitpc[aniter]);
+
+
+              if (ichg == 2) {
+                ckstartpc[aniter] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),nbinssh,0.0,shqmax);
+                anetaphitpc[aniter]->AddCorrFctn(ckstartpc[aniter]);
+              }
+              else {
+
+                cqinvtpc[aniter] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2*nbinssh,0.0,2*shqmax);
+                anetaphitpc[aniter]->AddCorrFctn(cqinvtpc[aniter]);
+
+              }
+
+              cylmtpc[aniter] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%i", chrgs[ichg], imult),2,nbinssh, 0.0,shqmax,runshlcms);
+              anetaphitpc[aniter]->AddCorrFctn(cylmtpc[aniter]);
+
+
+              // cAvgSeptpc[aniter] = new AliFemtoAvgSepCorrFctn(Form("cAvgSep%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),4*nbinssh,0.0,200);
+              // anetaphitpc[aniter]->AddCorrFctn(cAvgSeptpc[aniter]);
+
+              // cqinvinnertpc[aniter] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%d", chrgs[ichg], imult),nbinssh,0.0,shqmax);
+              // cqinvinnertpc[aniter]->SetRadius(1.2);
+              // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[aniter]);
+
+
+              if (runktdep) {
+                int ktm;
+                for (int ikt=0; ikt<numOfkTbins; ikt++) {
+
+                  ktm = aniter * numOfkTbins + ikt;
+                  ktpcuts[ktm] = new AliFemtoKTPairCut(ktrng[ikt], ktrng[ikt+1]);
+
+
+                  cylmkttpc[ktm] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%ikT%i", chrgs[ichg], imult, ikt),2,nbinssh,0.0,shqmax,runshlcms);
+                  cylmkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  anetaphitpc[aniter]->AddCorrFctn(cylmkttpc[ktm]);
+
+                  if (ichg == 2) {
+                    ckstarkttpc[ktm] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),nbinssh,0.0,shqmax);
+                    ckstarkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(ckstarkttpc[ktm]);
+
+
+                  }
+                  else {
+                    cqinvkttpc[ktm] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),2*nbinssh,0.0,2*shqmax);
+                    cqinvkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(cqinvkttpc[ktm]);
+                  }
+
+                  //         cqinvsqtpc[ktm] = new AliFemtoShareQualityCorrFctn(Form("cqinvsq%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  //         cqinvsqtpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //         anetaphitpc[aniter]->AddCorrFctn(cqinvsqtpc[ktm]);
+
+                  // cqinvinnertpc[ktm] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  // cqinvinnertpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  // cqinvinnertpc[ktm]->SetRadius(1.2);
+                  // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[ktm]);
+
+                  //         if (run3d) {
+                  //           cq3dlcmskttpc[ktm] = new AliFemtoCorrFctn3DLCMSSym(Form("cq3d%stpcM%ikT%i", chrgs[ichg], imult, ikt),60,(imult>3)?((imult>6)?((imult>7)?0.6:0.4):0.25):0.15);
+                  //           cq3dlcmskttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //           anetaphitpc[aniter]->AddCorrFctn(cq3dlcmskttpc[ktm]);
+                  //         }
+                }
+              }
+
+              // cdedpetaphi[aniter] = new AliFemtoCorrFctnDEtaDPhi(Form("cdedp%stpcM%i", chrgs[ichg], imult),240, 240);
+              // anetaphitpc[aniter]->AddCorrFctn(cdedpetaphi[aniter]);
+
+              Manager->AddAnalysis(anetaphitpc[aniter]);
+            }
+          }
+        }
+      }
+    }
+  }
+
+
+// *** End pion-pion (positive) analysis
+
+
+  return Manager;
+}
diff --git a/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcTof/monitors/ConfigFemtoAnalysis.C b/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcTof/monitors/ConfigFemtoAnalysis.C
new file mode 100644 (file)
index 0000000..1764e3e
--- /dev/null
@@ -0,0 +1,447 @@
+/*********************************************************************
+ *                                                                   *
+ * Configfemtoanalysis.C - configuration macro for the femtoscopic   *
+ * analysis, meant as a QA process for two-particle effects          *
+ *                                                                   *
+ * Author: Adam Kisiel (Adam.Kisiel@cern.ch)                         *
+ *                                                                   *
+ *********************************************************************/
+
+#if !defined(__CINT__) || defined(__MAKECINT_)
+#include "AliFemtoManager.h"
+#include "AliFemtoEventReaderESDChain.h"
+#include "AliFemtoEventReaderESDChainKine.h"
+#include "AliFemtoEventReaderAODChain.h"
+#include "AliFemtoSimpleAnalysis.h"
+#include "AliFemtoBasicEventCut.h"
+#include "AliFemtoESDTrackCut.h"
+#include "AliFemtoCorrFctn.h"
+#include "AliFemtoCutMonitorParticleYPt.h"
+#include "AliFemtoCutMonitorParticleVertPos.h"
+#include "AliFemtoCutMonitorParticleMomRes.h"
+#include "AliFemtoCutMonitorParticlePID.h"
+#include "AliFemtoCutMonitorEventMult.h"
+#include "AliFemtoCutMonitorEventVertex.h"
+#include "AliFemtoShareQualityTPCEntranceSepPairCut.h"
+#include "AliFemtoPairCutAntiGamma.h"
+#include "AliFemtoPairCutRadialDistance.h"
+#include "AliFemtoQinvCorrFctn.h"
+#include "AliFemtoCorrFctnNonIdDR.h"
+#include "AliFemtoShareQualityCorrFctn.h"
+#include "AliFemtoTPCInnerCorrFctn.h"
+#include "AliFemtoVertexMultAnalysis.h"
+#include "AliFemtoCorrFctn3DSpherical.h"
+#include "AliFemtoChi2CorrFctn.h"
+#include "AliFemtoCorrFctnTPCNcls.h"
+#include "AliFemtoBPLCMS3DCorrFctn.h"
+#include "AliFemtoCorrFctn3DLCMSSym.h"
+#include "AliFemtoModelBPLCMSCorrFctn.h"
+#include "AliFemtoModelCorrFctn3DSpherical.h"
+#include "AliFemtoModelGausLCMSFreezeOutGenerator.h"
+#include "AliFemtoModelGausRinvFreezeOutGenerator.h"
+#include "AliFemtoModelManager.h"
+#include "AliFemtoModelWeightGeneratorBasic.h"
+#include "AliFemtoModelWeightGeneratorLednicky.h"
+#include "AliFemtoCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnSource.h"
+#include "AliFemtoCutMonitorParticlePtPDG.h"
+#include "AliFemtoKTPairCut.h"
+#include "AliFemtoAvgSepCorrFctn.h"
+#endif
+
+//________________________________________________________________________
+AliFemtoManager* ConfigFemtoAnalysis() {
+
+  double PionMass = 0.13956995;
+  double KaonMass = 0.493677;
+  double ProtonMass = 0.938272013;
+
+  // double psi = TMath::Pi()/2.;
+  // double psid = TMath::Pi()/6.;
+
+  // int runepvzero[7] = {1, 1, 1, 1, 1, 1, 1};
+  // double epvzerobins[7] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid, -psi+4*psid, -psi+5*psid, -psi+6*psid};
+
+  double psi = TMath::Pi()/2.;
+  double psid = TMath::Pi()/3.;
+
+  int runepvzero[4] = {0, 0, 0, 1};
+  double epvzerobins[4] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid};
+
+  int runmults[10] = {1, 1, 0, 0, 0, 0, 0, 0, 0, 0};
+  int multbins[11] = {0.001, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900};
+
+  int runch[3] = {1, 1, 1};
+  const char *chrgs[3] = { "PP", "APAP", "PAP" };
+
+  int runktdep = 0;
+  double ktrng[3] = {0.01, 1.0, 5.0};
+
+  int numOfMultBins = 10;
+  int numOfChTypes = 3;
+  int numOfkTbins = 2;
+  int numOfEPvzero = 4;
+
+  int runqinv = 1;
+  int runshlcms = 0;// 0:PRF(PAP), 1:LCMS(PP,APAP)
+
+  int runtype = 2; // Types 0 - global, 1 - ITS only, 2 - TPC Inner
+  int isrealdata = 1;
+
+  //  int gammacut = 1;
+
+  double shqmax = 1.0;
+  int nbinssh = 100;
+
+  AliFemtoEventReaderAODChain *Reader = new AliFemtoEventReaderAODChain();
+  Reader->SetFilterBit(7);
+  Reader->SetCentralityPreSelection(0.001, 510);
+  Reader->SetEPVZERO(kTRUE);
+
+  AliFemtoManager* Manager = new AliFemtoManager();
+  Manager->SetEventReader(Reader);
+
+  AliFemtoVertexMultAnalysis    *anetaphitpc[10*3*2];
+  AliFemtoBasicEventCut         *mecetaphitpc[10*3*2];
+  AliFemtoCutMonitorEventMult   *cutPassEvMetaphitpc[50];
+  AliFemtoCutMonitorEventMult   *cutFailEvMetaphitpc[50];
+  AliFemtoCutMonitorEventVertex *cutPassEvVetaphitpc[50];
+  AliFemtoCutMonitorEventVertex *cutFailEvVetaphitpc[50];
+  AliFemtoESDTrackCut           *dtc1etaphitpc[50];
+  AliFemtoESDTrackCut           *dtc2etaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass2PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail2PIDetaphitpc[50];
+  //   AliFemtoPairCutAntiGamma      *sqpcetaphitpcdiff[10*3];
+  //   AliFemtoShareQualityTPCEntranceSepPairCut      *sqpcetaphitpcsame[10*3];
+  //AliFemtoPairCutAntiGamma      *sqpcetaphitpc[10*3];
+  AliFemtoPairCutRadialDistance      *sqpcetaphitpc[50];
+  //  AliFemtoChi2CorrFctn          *cchiqinvetaphitpc[20*2];
+  AliFemtoKTPairCut             *ktpcuts[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmtpc[50];
+  AliFemtoCorrFctnDirectYlm     *cylmkttpc[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmetaphitpc[10*3];
+  AliFemtoQinvCorrFctn          *cqinvkttpc[50*2];
+  AliFemtoQinvCorrFctn          *cqinvtpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstartpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstarkttpc[50*2];
+  AliFemtoCorrFctnDEtaDPhi      *cdedpetaphi[50*2];
+  AliFemtoAvgSepCorrFctn          *cAvgSeptpc[50];
+
+  //   AliFemtoCorrFctn3DLCMSSym     *cq3dlcmskttpc[20*2];
+  //   AliFemtoCorrFctnTPCNcls       *cqinvnclstpc[20];
+  //   AliFemtoShareQualityCorrFctn  *cqinvsqtpc[20*10];
+  //   AliFemtoChi2CorrFctn          *cqinvchi2tpc[20];
+    AliFemtoTPCInnerCorrFctn      *cqinvinnertpc[50];
+
+  // *** Third QA task - HBT analysis with all pair cuts off, TPC only ***
+  // *** Begin pion-pion (positive) analysis ***
+  int aniter = 0;
+
+  for (int imult = 0; imult < numOfMultBins; imult++) {
+    if (runmults[imult]) {
+
+      for (int ichg = 0; ichg < numOfChTypes; ichg++) {
+        if (runch[ichg]) {
+
+          for (int iepvzero = 0; iepvzero < numOfEPvzero; iepvzero++) {
+            if (runepvzero[iepvzero]) {
+
+              aniter = imult * numOfChTypes + ichg * numOfEPvzero + iepvzero;
+              // aniter = ichg * numOfMultBins + imult * numOfEPvzero + iepvzero;
+
+              // cout << "aniter = " << aniter << endl;
+              //        aniter = ichg * numOfMultBins + imult;
+
+              // if (ichg == 2)
+              //       runshlcms = 0;
+              // else
+              //       runshlcms = 1;
+
+
+              //________________________
+
+              anetaphitpc[aniter] = new AliFemtoVertexMultAnalysis(8, -8.0, 8.0, 4, multbins[imult], multbins[imult+1]);
+              anetaphitpc[aniter]->SetNumEventsToMix(10);
+              anetaphitpc[aniter]->SetMinSizePartCollection(1);
+              anetaphitpc[aniter]->SetVerboseMode(kFALSE);
+
+              mecetaphitpc[aniter] = new AliFemtoBasicEventCut();
+              mecetaphitpc[aniter]->SetEventMult(0.001,100000);
+              mecetaphitpc[aniter]->SetVertZPos(-8,8);
+
+              if (iepvzero == 3)
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[0],epvzerobins[3]);
+              else
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[iepvzero],epvzerobins[iepvzero+1]);
+
+              // if (isrealdata)
+              //     mecetaphitpc[aniter]->SetAcceptOnlyPhysics(kTRUE);
+
+              cutPassEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutPass%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              cutFailEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutFail%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              mecetaphitpc[aniter]->AddCutMonitor(cutPassEvMetaphitpc[aniter], cutFailEvMetaphitpc[aniter]);
+
+              cutPassEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutPass%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              cutFailEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutFail%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              mecetaphitpc[aniter]->AddCutMonitor(cutPassEvVetaphitpc[aniter], cutFailEvVetaphitpc[aniter]);
+
+              dtc1etaphitpc[aniter] = new AliFemtoESDTrackCut();
+              dtc2etaphitpc[aniter] = new AliFemtoESDTrackCut();
+
+              if (ichg == 0) {
+                dtc1etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 1) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 2) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc2etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+                dtc2etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+
+              dtc1etaphitpc[aniter]->SetEta(-0.8,0.8);
+              dtc1etaphitpc[aniter]->SetMass(ProtonMass);
+              dtc1etaphitpc[aniter]->SetMostProbableProton();
+              dtc1etaphitpc[aniter]->SetNsigma(3.0);
+              //dtc1etaphitpc[aniter]->SetNsigma(2.0);
+              dtc1etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+              //dtc1etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              if (ichg == 2) {
+                dtc2etaphitpc[aniter]->SetEta(-0.8,0.8);
+                dtc2etaphitpc[aniter]->SetMass(ProtonMass);
+                dtc2etaphitpc[aniter]->SetMostProbableProton();
+                dtc2etaphitpc[aniter]->SetNsigma(3.0);
+                //dtc2etaphitpc[aniter]->SetNsigma(2.0);
+                dtc2etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+                //dtc2etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              }
+
+              // Track quality cuts
+
+              if (runtype == 0) {
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                //            dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0182, 0.0350, -1.01);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 1) {
+                //      dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit|AliESDtrack::kITSpureSA);
+                //      dtc1etaphitpc[aniter]->SetminTPCncls(70);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                //      dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.25);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 2) {
+                //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                //dtc1etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                // dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+
+                if (ichg == 2) {
+                  //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                  dtc2etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                  dtc2etaphitpc[aniter]->SetminTPCncls(80);
+                  dtc2etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                  dtc2etaphitpc[aniter]->SetLabel(kFALSE);
+                  dtc2etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                  //dtc2etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                  // dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+                }
+
+              }
+
+              cutPass1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              cutFail1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1YPtetaphitpc[aniter], cutFail1YPtetaphitpc[aniter]);
+
+              cutPass1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2);//0-pion,1-kaon,2-proton
+              cutFail1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult , iepvzero),2);
+              dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1PIDetaphitpc[aniter], cutFail1PIDetaphitpc[aniter]);
+
+              // if (ichg == 2){
+              //     cutPass2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass2%stpcM%i", chrgs[ichg], imult),2);//0-pion,1-kaon,2-proton
+              //     cutFail2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail2%stpcM%i", chrgs[ichg], imult),2);
+              //     dtc2etaphitpc[aniter]->AddCutMonitor(cutPass2PIDetaphitpc[aniter], cutFail2PIDetaphitpc[aniter]);
+              // }
+
+              // sqpcetaphitpc[aniter] = new AliFemtoPairCutAntiGamma();
+              sqpcetaphitpc[aniter] = new AliFemtoPairCutRadialDistance();
+
+              if (runtype == 0) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(1.5);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(0.12, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 1) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(1.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(5.0);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(1.2, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 2) {
+                //sqpcetaphitpc[aniter]->SetUseAOD(kTRUE);
+
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+
+                //         if (gammacut == 0) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //}
+                //else if (gammacut == 1) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //}
+
+                // sqpcetaphitpc[aniter]->SetMagneticFieldSign(-1); // field1 -1, field3 +1
+                // sqpcetaphitpc[aniter]->SetMinimumRadius(0.8); // biggest inefficiency for R=1.1 m (checked on small sample)
+
+                sqpcetaphitpc[aniter]->SetMinimumRadius(1.2); //0.8
+                sqpcetaphitpc[aniter]->SetPhiStarMin(kFALSE);
+                sqpcetaphitpc[aniter]->SetPhiStarDifferenceMinimum(0.045); // 0.012 - pions, 0.017 - kaons, 0.018
+                sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.01); // 0.017 - pions, 0.015 - kaons
+
+              }
+
+              anetaphitpc[aniter]->SetEventCut(mecetaphitpc[aniter]);
+
+              if (ichg == 2) {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc2etaphitpc[aniter]);
+              }
+              else {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc1etaphitpc[aniter]);
+              }
+
+              anetaphitpc[aniter]->SetPairCut(sqpcetaphitpc[aniter]);
+
+
+              // if (ichg == 2) {
+              //   ckstartpc[aniter] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),nbinssh,0.0,shqmax);
+              //   anetaphitpc[aniter]->AddCorrFctn(ckstartpc[aniter]);
+              // }
+              // else {
+
+              //   cqinvtpc[aniter] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2*nbinssh,0.0,2*shqmax);
+              //   anetaphitpc[aniter]->AddCorrFctn(cqinvtpc[aniter]);
+
+              // }
+
+              // cylmtpc[aniter] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%i", chrgs[ichg], imult),2,nbinssh, 0.0,shqmax,runshlcms);
+              // anetaphitpc[aniter]->AddCorrFctn(cylmtpc[aniter]);
+
+
+              // cAvgSeptpc[aniter] = new AliFemtoAvgSepCorrFctn(Form("cAvgSep%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),4*nbinssh,0.0,200);
+              // anetaphitpc[aniter]->AddCorrFctn(cAvgSeptpc[aniter]);
+
+              // cqinvinnertpc[aniter] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%d", chrgs[ichg], imult),nbinssh,0.0,shqmax);
+              // cqinvinnertpc[aniter]->SetRadius(1.2);
+              // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[aniter]);
+
+
+              if (runktdep) {
+                int ktm;
+                for (int ikt=0; ikt<numOfkTbins; ikt++) {
+
+                  ktm = aniter * numOfkTbins + ikt;
+                  ktpcuts[ktm] = new AliFemtoKTPairCut(ktrng[ikt], ktrng[ikt+1]);
+
+
+                  cylmkttpc[ktm] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%ikT%i", chrgs[ichg], imult, ikt),2,nbinssh,0.0,shqmax,runshlcms);
+                  cylmkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  anetaphitpc[aniter]->AddCorrFctn(cylmkttpc[ktm]);
+
+                  if (ichg == 2) {
+                    ckstarkttpc[ktm] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),nbinssh,0.0,shqmax);
+                    ckstarkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(ckstarkttpc[ktm]);
+
+
+                  }
+                  else {
+                    cqinvkttpc[ktm] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),2*nbinssh,0.0,2*shqmax);
+                    cqinvkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(cqinvkttpc[ktm]);
+                  }
+
+                  //         cqinvsqtpc[ktm] = new AliFemtoShareQualityCorrFctn(Form("cqinvsq%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  //         cqinvsqtpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //         anetaphitpc[aniter]->AddCorrFctn(cqinvsqtpc[ktm]);
+
+                  // cqinvinnertpc[ktm] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  // cqinvinnertpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  // cqinvinnertpc[ktm]->SetRadius(1.2);
+                  // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[ktm]);
+
+                  //         if (run3d) {
+                  //           cq3dlcmskttpc[ktm] = new AliFemtoCorrFctn3DLCMSSym(Form("cq3d%stpcM%ikT%i", chrgs[ichg], imult, ikt),60,(imult>3)?((imult>6)?((imult>7)?0.6:0.4):0.25):0.15);
+                  //           cq3dlcmskttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //           anetaphitpc[aniter]->AddCorrFctn(cq3dlcmskttpc[ktm]);
+                  //         }
+                }
+              }
+
+              // cdedpetaphi[aniter] = new AliFemtoCorrFctnDEtaDPhi(Form("cdedp%stpcM%i", chrgs[ichg], imult),240, 240);
+              // anetaphitpc[aniter]->AddCorrFctn(cdedpetaphi[aniter]);
+
+              Manager->AddAnalysis(anetaphitpc[aniter]);
+            }
+          }
+        }
+      }
+    }
+  }
+
+
+// *** End pion-pion (positive) analysis
+
+
+  return Manager;
+}
diff --git a/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcTof/semicentral/ConfigFemtoAnalysis.C b/PWGCF/FEMTOSCOPY/macros/Train/ProtonFemtoscopy/PidTpcTof/semicentral/ConfigFemtoAnalysis.C
new file mode 100644 (file)
index 0000000..82fdb9f
--- /dev/null
@@ -0,0 +1,447 @@
+/*********************************************************************
+ *                                                                   *
+ * Configfemtoanalysis.C - configuration macro for the femtoscopic   *
+ * analysis, meant as a QA process for two-particle effects          *
+ *                                                                   *
+ * Author: Adam Kisiel (Adam.Kisiel@cern.ch)                         *
+ *                                                                   *
+ *********************************************************************/
+
+#if !defined(__CINT__) || defined(__MAKECINT_)
+#include "AliFemtoManager.h"
+#include "AliFemtoEventReaderESDChain.h"
+#include "AliFemtoEventReaderESDChainKine.h"
+#include "AliFemtoEventReaderAODChain.h"
+#include "AliFemtoSimpleAnalysis.h"
+#include "AliFemtoBasicEventCut.h"
+#include "AliFemtoESDTrackCut.h"
+#include "AliFemtoCorrFctn.h"
+#include "AliFemtoCutMonitorParticleYPt.h"
+#include "AliFemtoCutMonitorParticleVertPos.h"
+#include "AliFemtoCutMonitorParticleMomRes.h"
+#include "AliFemtoCutMonitorParticlePID.h"
+#include "AliFemtoCutMonitorEventMult.h"
+#include "AliFemtoCutMonitorEventVertex.h"
+#include "AliFemtoShareQualityTPCEntranceSepPairCut.h"
+#include "AliFemtoPairCutAntiGamma.h"
+#include "AliFemtoPairCutRadialDistance.h"
+#include "AliFemtoQinvCorrFctn.h"
+#include "AliFemtoCorrFctnNonIdDR.h"
+#include "AliFemtoShareQualityCorrFctn.h"
+#include "AliFemtoTPCInnerCorrFctn.h"
+#include "AliFemtoVertexMultAnalysis.h"
+#include "AliFemtoCorrFctn3DSpherical.h"
+#include "AliFemtoChi2CorrFctn.h"
+#include "AliFemtoCorrFctnTPCNcls.h"
+#include "AliFemtoBPLCMS3DCorrFctn.h"
+#include "AliFemtoCorrFctn3DLCMSSym.h"
+#include "AliFemtoModelBPLCMSCorrFctn.h"
+#include "AliFemtoModelCorrFctn3DSpherical.h"
+#include "AliFemtoModelGausLCMSFreezeOutGenerator.h"
+#include "AliFemtoModelGausRinvFreezeOutGenerator.h"
+#include "AliFemtoModelManager.h"
+#include "AliFemtoModelWeightGeneratorBasic.h"
+#include "AliFemtoModelWeightGeneratorLednicky.h"
+#include "AliFemtoCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnDirectYlm.h"
+#include "AliFemtoModelCorrFctnSource.h"
+#include "AliFemtoCutMonitorParticlePtPDG.h"
+#include "AliFemtoKTPairCut.h"
+#include "AliFemtoAvgSepCorrFctn.h"
+#endif
+
+//________________________________________________________________________
+AliFemtoManager* ConfigFemtoAnalysis() {
+
+  double PionMass = 0.13956995;
+  double KaonMass = 0.493677;
+  double ProtonMass = 0.938272013;
+
+  // double psi = TMath::Pi()/2.;
+  // double psid = TMath::Pi()/6.;
+
+  // int runepvzero[7] = {1, 1, 1, 1, 1, 1, 1};
+  // double epvzerobins[7] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid, -psi+4*psid, -psi+5*psid, -psi+6*psid};
+
+  double psi = TMath::Pi()/2.;
+  double psid = TMath::Pi()/3.;
+
+  int runepvzero[4] = {0, 0, 0, 1};
+  double epvzerobins[4] = {-psi, -psi+psid, -psi+2*psid, -psi+3*psid};
+
+  int runmults[10] = {0, 0, 1, 1, 1, 1, 0, 0, 0, 0};
+  int multbins[11] = {0.001, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900};
+
+  int runch[3] = {1, 1, 1};
+  const char *chrgs[3] = { "PP", "APAP", "PAP" };
+
+  int runktdep = 1;
+  double ktrng[3] = {0.01, 1.0, 5.0};
+
+  int numOfMultBins = 10;
+  int numOfChTypes = 3;
+  int numOfkTbins = 2;
+  int numOfEPvzero = 4;
+
+  int runqinv = 1;
+  int runshlcms = 0;// 0:PRF(PAP), 1:LCMS(PP,APAP)
+
+  int runtype = 2; // Types 0 - global, 1 - ITS only, 2 - TPC Inner
+  int isrealdata = 1;
+
+  //  int gammacut = 1;
+
+  double shqmax = 1.0;
+  int nbinssh = 100;
+
+  AliFemtoEventReaderAODChain *Reader = new AliFemtoEventReaderAODChain();
+  Reader->SetFilterBit(7);
+  Reader->SetCentralityPreSelection(0.001, 510);
+  Reader->SetEPVZERO(kTRUE);
+
+  AliFemtoManager* Manager = new AliFemtoManager();
+  Manager->SetEventReader(Reader);
+
+  AliFemtoVertexMultAnalysis    *anetaphitpc[10*3*2];
+  AliFemtoBasicEventCut         *mecetaphitpc[10*3*2];
+  AliFemtoCutMonitorEventMult   *cutPassEvMetaphitpc[50];
+  AliFemtoCutMonitorEventMult   *cutFailEvMetaphitpc[50];
+  // AliFemtoCutMonitorEventVertex *cutPassEvVetaphitpc[50];
+  // AliFemtoCutMonitorEventVertex *cutFailEvVetaphitpc[50];
+  AliFemtoESDTrackCut           *dtc1etaphitpc[50];
+  AliFemtoESDTrackCut           *dtc2etaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail1YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail1PIDetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutPass2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticleYPt *cutFail2YPtetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutPass2PIDetaphitpc[50];
+  AliFemtoCutMonitorParticlePID *cutFail2PIDetaphitpc[50];
+  //   AliFemtoPairCutAntiGamma      *sqpcetaphitpcdiff[10*3];
+  //   AliFemtoShareQualityTPCEntranceSepPairCut      *sqpcetaphitpcsame[10*3];
+  //AliFemtoPairCutAntiGamma      *sqpcetaphitpc[10*3];
+  AliFemtoPairCutRadialDistance      *sqpcetaphitpc[50];
+  //  AliFemtoChi2CorrFctn          *cchiqinvetaphitpc[20*2];
+  AliFemtoKTPairCut             *ktpcuts[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmtpc[50];
+  AliFemtoCorrFctnDirectYlm     *cylmkttpc[50*2];
+  AliFemtoCorrFctnDirectYlm     *cylmetaphitpc[10*3];
+  AliFemtoQinvCorrFctn          *cqinvkttpc[50*2];
+  AliFemtoQinvCorrFctn          *cqinvtpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstartpc[50];
+  AliFemtoCorrFctnNonIdDR       *ckstarkttpc[50*2];
+  AliFemtoCorrFctnDEtaDPhi      *cdedpetaphi[50*2];
+  AliFemtoAvgSepCorrFctn          *cAvgSeptpc[50];
+
+  //   AliFemtoCorrFctn3DLCMSSym     *cq3dlcmskttpc[20*2];
+  //   AliFemtoCorrFctnTPCNcls       *cqinvnclstpc[20];
+  //   AliFemtoShareQualityCorrFctn  *cqinvsqtpc[20*10];
+  //   AliFemtoChi2CorrFctn          *cqinvchi2tpc[20];
+    AliFemtoTPCInnerCorrFctn      *cqinvinnertpc[50];
+
+  // *** Third QA task - HBT analysis with all pair cuts off, TPC only ***
+  // *** Begin pion-pion (positive) analysis ***
+  int aniter = 0;
+
+  for (int imult = 0; imult < numOfMultBins; imult++) {
+    if (runmults[imult]) {
+
+      for (int ichg = 0; ichg < numOfChTypes; ichg++) {
+        if (runch[ichg]) {
+
+          for (int iepvzero = 0; iepvzero < numOfEPvzero; iepvzero++) {
+            if (runepvzero[iepvzero]) {
+
+              aniter = imult * numOfChTypes + ichg * numOfEPvzero + iepvzero;
+              // aniter = ichg * numOfMultBins + imult * numOfEPvzero + iepvzero;
+
+              // cout << "aniter = " << aniter << endl;
+              //        aniter = ichg * numOfMultBins + imult;
+
+              // if (ichg == 2)
+              //       runshlcms = 0;
+              // else
+              //       runshlcms = 1;
+
+
+              //________________________
+
+              anetaphitpc[aniter] = new AliFemtoVertexMultAnalysis(8, -8.0, 8.0, 4, multbins[imult], multbins[imult+1]);
+              anetaphitpc[aniter]->SetNumEventsToMix(10);
+              anetaphitpc[aniter]->SetMinSizePartCollection(1);
+              anetaphitpc[aniter]->SetVerboseMode(kFALSE);
+
+              mecetaphitpc[aniter] = new AliFemtoBasicEventCut();
+              mecetaphitpc[aniter]->SetEventMult(0.001,100000);
+              mecetaphitpc[aniter]->SetVertZPos(-8,8);
+
+              if (iepvzero == 3)
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[0],epvzerobins[3]);
+              else
+                mecetaphitpc[aniter]->SetEPVZERO(epvzerobins[iepvzero],epvzerobins[iepvzero+1]);
+
+              // if (isrealdata)
+              //     mecetaphitpc[aniter]->SetAcceptOnlyPhysics(kTRUE);
+
+              // cutPassEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutPass%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              // cutFailEvMetaphitpc[aniter] = new AliFemtoCutMonitorEventMult(Form("cutFail%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero));
+              // mecetaphitpc[aniter]->AddCutMonitor(cutPassEvMetaphitpc[aniter], cutFailEvMetaphitpc[aniter]);
+
+              // cutPassEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutPass%stpcM%i", chrgs[ichg], imult));
+              // cutFailEvVetaphitpc[aniter] = new AliFemtoCutMonitorEventVertex(Form("cutFail%stpcM%i", chrgs[ichg], imult));
+              // mecetaphitpc[aniter]->AddCutMonitor(cutPassEvVetaphitpc[aniter], cutFailEvVetaphitpc[aniter]);
+
+              dtc1etaphitpc[aniter] = new AliFemtoESDTrackCut();
+              dtc2etaphitpc[aniter] = new AliFemtoESDTrackCut();
+
+              if (ichg == 0) {
+                dtc1etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 1) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+              else if (ichg == 2) {
+                dtc1etaphitpc[aniter]->SetCharge(-1.0);
+                dtc2etaphitpc[aniter]->SetCharge(1.0);
+                dtc1etaphitpc[aniter]->SetPt(0.7,4.0);
+                dtc2etaphitpc[aniter]->SetPt(0.7,4.0);
+              }
+
+              dtc1etaphitpc[aniter]->SetEta(-0.8,0.8);
+              dtc1etaphitpc[aniter]->SetMass(ProtonMass);
+              dtc1etaphitpc[aniter]->SetMostProbableProton();
+              dtc1etaphitpc[aniter]->SetNsigma(3.0);
+              //dtc1etaphitpc[aniter]->SetNsigma(2.0);
+              dtc1etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+              //dtc1etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              if (ichg == 2) {
+                dtc2etaphitpc[aniter]->SetEta(-0.8,0.8);
+                dtc2etaphitpc[aniter]->SetMass(ProtonMass);
+                dtc2etaphitpc[aniter]->SetMostProbableProton();
+                dtc2etaphitpc[aniter]->SetNsigma(3.0);
+                //dtc2etaphitpc[aniter]->SetNsigma(2.0);
+                dtc2etaphitpc[aniter]->SetNsigmaTPCTOF(kTRUE);
+                //dtc2etaphitpc[aniter]->SetNsigmaTPConly(kTRUE);
+
+              }
+
+              // Track quality cuts
+
+              if (runtype == 0) {
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                //            dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0182, 0.0350, -1.01);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 1) {
+                //      dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                //    dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit);
+                //         dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit|AliESDtrack::kITSpureSA);
+                //      dtc1etaphitpc[aniter]->SetminTPCncls(70);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                //    dtc1etaphitpc[aniter]->SetMaxITSChiNdof(6.0);
+                //      dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(6.0);
+                dtc1etaphitpc[aniter]->SetMaxImpactXY(0.2);
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.25);
+                //      dtc1etaphitpc[aniter]->SetMaxSigmaToVertex(6.0);
+              }
+              else if (runtype == 2) {
+                //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                dtc1etaphitpc[aniter]->SetminTPCncls(80);
+                dtc1etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                dtc1etaphitpc[aniter]->SetLabel(kFALSE);
+                dtc1etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                //dtc1etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                // dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                dtc1etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+
+                if (ichg == 2) {
+                  //dtc1etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
+                  dtc2etaphitpc[aniter]->SetStatus(AliESDtrack::kTPCin);
+                  dtc2etaphitpc[aniter]->SetminTPCncls(80);
+                  dtc2etaphitpc[aniter]->SetRemoveKinks(kTRUE);
+                  dtc2etaphitpc[aniter]->SetLabel(kFALSE);
+                  dtc2etaphitpc[aniter]->SetMaxTPCChiNdof(4.0);
+                  //dtc2etaphitpc[aniter]->SetMaxImpactXY(0.1); // 2.4 0.1
+                  // dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.0205, 0.035, -1.1);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactXYPtDep(0.018, 0.035, -1.01);     //      DCA xy
+                  dtc2etaphitpc[aniter]->SetMaxImpactZ(0.15); // 2.0 0.1
+                }
+
+              }
+
+              // cutPass1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              // cutFail1YPtetaphitpc[aniter] = new AliFemtoCutMonitorParticleYPt(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),ProtonMass);
+              // dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1YPtetaphitpc[aniter], cutFail1YPtetaphitpc[aniter]);
+
+              // cutPass1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass1%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2);//0-pion,1-kaon,2-proton
+              // cutFail1PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail1%stpcM%iPsi%i", chrgs[ichg], imult , iepvzero),2);
+              // dtc1etaphitpc[aniter]->AddCutMonitor(cutPass1PIDetaphitpc[aniter], cutFail1PIDetaphitpc[aniter]);
+
+              // if (ichg == 2){
+              //     cutPass2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutPass2%stpcM%i", chrgs[ichg], imult),2);//0-pion,1-kaon,2-proton
+              //     cutFail2PIDetaphitpc[aniter] = new AliFemtoCutMonitorParticlePID(Form("cutFail2%stpcM%i", chrgs[ichg], imult),2);
+              //     dtc2etaphitpc[aniter]->AddCutMonitor(cutPass2PIDetaphitpc[aniter], cutFail2PIDetaphitpc[aniter]);
+              // }
+
+              // sqpcetaphitpc[aniter] = new AliFemtoPairCutAntiGamma();
+              sqpcetaphitpc[aniter] = new AliFemtoPairCutRadialDistance();
+
+              if (runtype == 0) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(1.5);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(0.12, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 1) {
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(1.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+                // sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                // sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //         sqpcetaphitpc[aniter]->SetTPCEntranceSepMinimum(5.0);
+                //sqpcetaphitpc[aniter]->SetRadialDistanceMinimum(1.2, 0.03);
+                //         sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.02);
+              }
+              else if (runtype == 2) {
+                //sqpcetaphitpc[aniter]->SetUseAOD(kTRUE);
+
+                sqpcetaphitpc[aniter]->SetShareQualityMax(1.0);
+                sqpcetaphitpc[aniter]->SetShareFractionMax(0.05);
+                sqpcetaphitpc[aniter]->SetRemoveSameLabel(kFALSE);
+
+                //         if (gammacut == 0) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.0);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.0);
+                //}
+                //else if (gammacut == 1) {
+                //sqpcetaphitpc[aniter]->SetMaxEEMinv(0.002);
+                //sqpcetaphitpc[aniter]->SetMaxThetaDiff(0.008);
+                //}
+
+                // sqpcetaphitpc[aniter]->SetMagneticFieldSign(-1); // field1 -1, field3 +1
+                // sqpcetaphitpc[aniter]->SetMinimumRadius(0.8); // biggest inefficiency for R=1.1 m (checked on small sample)
+
+                sqpcetaphitpc[aniter]->SetMinimumRadius(1.2); //0.8
+                sqpcetaphitpc[aniter]->SetPhiStarMin(kFALSE);
+                sqpcetaphitpc[aniter]->SetPhiStarDifferenceMinimum(0.045); // 0.012 - pions, 0.017 - kaons, 0.018
+                sqpcetaphitpc[aniter]->SetEtaDifferenceMinimum(0.01); // 0.017 - pions, 0.015 - kaons
+
+              }
+
+              anetaphitpc[aniter]->SetEventCut(mecetaphitpc[aniter]);
+
+              if (ichg == 2) {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc2etaphitpc[aniter]);
+              }
+              else {
+                anetaphitpc[aniter]->SetFirstParticleCut(dtc1etaphitpc[aniter]);
+                anetaphitpc[aniter]->SetSecondParticleCut(dtc1etaphitpc[aniter]);
+              }
+
+              anetaphitpc[aniter]->SetPairCut(sqpcetaphitpc[aniter]);
+
+
+              if (ichg == 2) {
+                ckstartpc[aniter] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),nbinssh,0.0,shqmax);
+                anetaphitpc[aniter]->AddCorrFctn(ckstartpc[aniter]);
+              }
+              else {
+
+                cqinvtpc[aniter] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),2*nbinssh,0.0,2*shqmax);
+                anetaphitpc[aniter]->AddCorrFctn(cqinvtpc[aniter]);
+
+              }
+
+              cylmtpc[aniter] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%i", chrgs[ichg], imult),2,nbinssh, 0.0,shqmax,runshlcms);
+              anetaphitpc[aniter]->AddCorrFctn(cylmtpc[aniter]);
+
+
+              // cAvgSeptpc[aniter] = new AliFemtoAvgSepCorrFctn(Form("cAvgSep%stpcM%iPsi%i", chrgs[ichg], imult, iepvzero),4*nbinssh,0.0,200);
+              // anetaphitpc[aniter]->AddCorrFctn(cAvgSeptpc[aniter]);
+
+              // cqinvinnertpc[aniter] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%d", chrgs[ichg], imult),nbinssh,0.0,shqmax);
+              // cqinvinnertpc[aniter]->SetRadius(1.2);
+              // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[aniter]);
+
+
+              if (runktdep) {
+                int ktm;
+                for (int ikt=0; ikt<numOfkTbins; ikt++) {
+
+                  ktm = aniter * numOfkTbins + ikt;
+                  ktpcuts[ktm] = new AliFemtoKTPairCut(ktrng[ikt], ktrng[ikt+1]);
+
+
+                  cylmkttpc[ktm] = new AliFemtoCorrFctnDirectYlm(Form("cylm%stpcM%ikT%i", chrgs[ichg], imult, ikt),2,nbinssh,0.0,shqmax,runshlcms);
+                  cylmkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  anetaphitpc[aniter]->AddCorrFctn(cylmkttpc[ktm]);
+
+                  if (ichg == 2) {
+                    ckstarkttpc[ktm] = new AliFemtoCorrFctnNonIdDR(Form("ckstar%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),nbinssh,0.0,shqmax);
+                    ckstarkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(ckstarkttpc[ktm]);
+
+
+                  }
+                  else {
+                    cqinvkttpc[ktm] = new AliFemtoQinvCorrFctn(Form("cqinv%stpcM%iPsi%ikT%i", chrgs[ichg], imult, iepvzero, ikt),2*nbinssh,0.0,2*shqmax);
+                    cqinvkttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                    anetaphitpc[aniter]->AddCorrFctn(cqinvkttpc[ktm]);
+                  }
+
+                  //         cqinvsqtpc[ktm] = new AliFemtoShareQualityCorrFctn(Form("cqinvsq%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  //         cqinvsqtpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //         anetaphitpc[aniter]->AddCorrFctn(cqinvsqtpc[ktm]);
+
+                  // cqinvinnertpc[ktm] = new AliFemtoTPCInnerCorrFctn(Form("cqinvinner%stpcM%ikT%i", chrgs[ichg], imult, ikt),nbinssh,0.0,shqmax);
+                  // cqinvinnertpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  // cqinvinnertpc[ktm]->SetRadius(1.2);
+                  // anetaphitpc[aniter]->AddCorrFctn(cqinvinnertpc[ktm]);
+
+                  //         if (run3d) {
+                  //           cq3dlcmskttpc[ktm] = new AliFemtoCorrFctn3DLCMSSym(Form("cq3d%stpcM%ikT%i", chrgs[ichg], imult, ikt),60,(imult>3)?((imult>6)?((imult>7)?0.6:0.4):0.25):0.15);
+                  //           cq3dlcmskttpc[ktm]->SetPairSelectionCut(ktpcuts[ktm]);
+                  //           anetaphitpc[aniter]->AddCorrFctn(cq3dlcmskttpc[ktm]);
+                  //         }
+                }
+              }
+
+              // cdedpetaphi[aniter] = new AliFemtoCorrFctnDEtaDPhi(Form("cdedp%stpcM%i", chrgs[ichg], imult),240, 240);
+              // anetaphitpc[aniter]->AddCorrFctn(cdedpetaphi[aniter]);
+
+              Manager->AddAnalysis(anetaphitpc[aniter]);
+            }
+          }
+        }
+      }
+    }
+  }
+
+
+// *** End pion-pion (positive) analysis
+
+
+  return Manager;
+}