TPC/AliESDcosmic.cxx

   1 #include "TMath.h"
   2 #include "TArrayI.h"
   3 #include "TDatabasePDG.h"
   4 #include "TTreeStream.h"
   5 #include "TVectorD.h"
   6
   7
   8 #include "AliESDEvent.h"
   9 #include "AliESDtrack.h"
  10 #include "AliTracker.h"
  11 #include "AliESDcosmic.h"
  12
  13 ClassImp(AliESDcosmic)
  14
  15 AliESDcosmic::AliESDcosmic():
  16     TObject(),
  17     fESD(0),
  18     fTracks(0),
  19     fTracksAcorde(0),
  20     fPair(0),
  21     fDebugStreamer(0)
  22 {
  23   //
  24   //
  25   //
  26 }
  27
  28
  29 AliESDcosmic::~AliESDcosmic(){
  30   //
  31   //
  32   //
  33   delete fPair;
  34   delete fTracksAcorde;
  35
  36 }
  37
  38 void AliESDcosmic::ProcessEvent(AliESDEvent* event){
  39   //
  40   //
  41   //
  42   fESD=event;
  43   TString  string =  fESD->GetFiredTriggerClasses();
  44   if (!(string.Contains("ASL"))) return;
  45   //
  46   const Float_t kCutMinDir=0.95;
  47   const Float_t kMaxD=1000.;
  48   Int_t ntracks=event->GetNumberOfTracks();
  49   if (ntracks<=0) return;
  50   if (!fTracks) fTracks= new TClonesArray("AliESDtrack",0);
  51   fPair = new TArrayI(ntracks);
  52
  53   for (Int_t i=0;i<ntracks;++i) {
  54     AliESDtrack *track0 = event->GetTrack(i);
  55     AliESDtrack *trackPair=0;
  56     (*fPair)[i]=-1;
  57     // track0 - choosen upper part
  58     if (!track0) continue;
  59     if (!track0->GetOuterParam()) continue;
  60     Double_t dir0[3];
  61     track0->GetDirection(dir0);
  62     Float_t minDist=kMaxD;
  63     //
  64     for (Int_t j=i;j<ntracks;++j) {
  65       AliESDtrack *track1 = event->GetTrack(j);
  66       //track 1 lower part
  67       if (!track1) continue;
  68       if (!track1->GetOuterParam()) continue;
  69       Double_t dir1[3];
  70       track1->GetDirection(dir1);
  71       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
  72       if (TMath::Abs(dir)<kCutMinDir) continue;               // direction vector product
  73       //
  74       // calculate distance
  75       Float_t dy = (track0->GetY()+track1->GetY());
  76       Float_t sy2 = track0->GetSigmaY2()+track1->GetSigmaY2();
  77       Float_t dphi = (track0->GetAlpha()-track1->GetAlpha()-TMath::Pi());
  78       Float_t sphi2  = track0->GetSigmaSnp2()+track1->GetSigmaSnp2();
  79       Float_t dtheta = (track0->GetTgl()-track1->GetTgl());
  80       Float_t stheta2 = track0->GetSigmaTgl2()+track1->GetSigmaTgl2();
  81       Float_t normDist = TMath::Sqrt(dy*dy/sy2+dphi*dphi/sphi2+dtheta*dtheta/stheta2);
  82       //
  83       if (normDist>minDist) continue;
  84       minDist = normDist;
  85       trackPair=track1;
  86       (*fPair)[i]=j;
  87     }
  88     //
  89   }
  90   PropagateToAcorde();
  91 }
  92
  93
  94 void  AliESDcosmic::PropagateToAcorde(){
  95   //
  96   //
  97   //
  98   const Double_t kRL3=510;   // radius of L3 magnet
  99   const Double_t kxAcorde=850.;
 100   //
 101   if (!fTracksAcorde) fTracksAcorde= new TClonesArray("AliExternalTrackParam",0);
 102   Int_t ntracks=fESD->GetNumberOfTracks();
 103   Int_t counter=0;
 104   //
 105   for (Int_t i=0; i<ntracks;i++){
 106     Int_t index = i;
 107     AliESDtrack * upperTrack = fESD->GetTrack(i);
 108     if (upperTrack->GetOuterParam()==0) continue;
 109     Double_t gxyz[3];
 110     upperTrack->GetOuterParam()->GetXYZ(gxyz);
 111     if ((*fPair)[i]>0){
 112       AliESDtrack * track2 = fESD->GetTrack((*fPair)[i]);
 113       if (track2->GetOuterParam()){
 114         Double_t gxyz2[3];
 115         track2->GetOuterParam()->GetXYZ(gxyz2);
 116         if (gxyz2[1]>gxyz[1]) {
 117           upperTrack=track2;
 118           index=(*fPair)[i];
 119         }
 120       }
 121     }
 122     //
 123     AliExternalTrackParam *upper = (AliExternalTrackParam *)(upperTrack->GetOuterParam()->Clone());
 124     Bool_t isOK = upper->PropagateTo(kRL3,fESD->GetMagneticField());
 125     upper->GetXYZ(gxyz);
 126     if (gxyz[1]<0) continue;
 127     for (Int_t iter=0; iter<20;iter++){
 128       upper->GetXYZ(gxyz);
 129       Double_t galpha = TMath::ATan2(gxyz[1],gxyz[0]);
 130       Double_t alpha=galpha;
 131       galpha*=180/TMath::Pi();
 132       if (iter>1){
 133         if (galpha<45.)   alpha = TMath::Pi()/8;
 134         if (galpha>135.) alpha =  TMath::Pi()*(1-1/8.);
 135         if (galpha>45.&&galpha<135.) alpha = TMath::Pi()/2.;
 136       }
 137       if (isOK) upper->Rotate(alpha);
 138       if (isOK) isOK = upper->PropagateTo(kxAcorde,0);
 139     }
 140     if (isOK) {
 141       new ((*fTracksAcorde)[index]) AliExternalTrackParam(*upper);
 142       counter++;
 143     }
 144   }
 145 }
 146
 147 void AliESDcosmic::DumpToTree(){
 148   //
 149   //
 150   //
 151   TTreeSRedirector * cstream = fDebugStreamer;
 152   if (!cstream) return;
 153   Int_t ntracks0 =fESD->GetNumberOfTracks();
 154   Int_t ntracks  = (fTracksAcorde) ? fTracksAcorde->GetEntries():0;
 155   Float_t mag = fESD->GetMagneticField();
 156   Int_t   run = fESD->GetRunNumber();
 157   Int_t   event = fESD->GetEventNumberInFile();
 158   AliESDHeader* header = fESD->GetHeader();
 159
 160   (*cstream)<<"eventInfo"<<
 161     "run="<<run<<
 162     "event="<<event<<
 163     "header="<<header<<
 164     "mag="<<mag<<
 165     "ntracks0="<<ntracks0<<
 166     "ntracks="<<ntracks<<
 167     "\n";
 168
 169   for (Int_t i=0;i<ntracks;i++){
 170     if (!fTracksAcorde->At(i)) continue;
 171     TVectorD gxyz(3);
 172     TVectorD gpxyz(3);
 173     AliExternalTrackParam * param = (AliExternalTrackParam *)fTracksAcorde->At(i);
 174     param->GetXYZ(gxyz.GetMatrixArray());
 175     param->GetPxPyPz(gpxyz.GetMatrixArray());
 176     (*cstream) << "esdCosmic" <<
 177       "run="<<run<<
 178       "event="<<event<<
 179       "header="<<header<<
 180       "mag="<<mag<<
 181       "ntracks0="<<ntracks0<<
 182       "ntracks="<<ntracks<<
 183       "tr.="<<param<<
 184       "p.="<<&gxyz<<
 185       "m.="<<&gpxyz<<
 186       "\n";
 187   }
 188
 189 }
 190