TRD/qaRec/AliTRDmultiplicity.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ////////////////////////////////////////////////////////////////////////////
  18 //                                                                        //
  19 //  TRD multiplicity                                                      //
  20 //                                                                        //
  21 // Task to select true single tracks
  22 //
  23 // Program flow:
  24 // Part A - task AliTRDmultiplicity (within the framework qaRec):
  25 //   For TRD standalone or TRD&TOF tracks I make a FitRiemanTilt
  26 //   The resulting parameterization is dumped into a DebugStreamer written to a file
  27 //
  28 // Part B – $ALICE_ROOT/TRD/qaRec/macros/TrackletsinTRD.C[h] (analysis macro):
  29 //   The clusters are read in
  30 //   The above produced file is read in
  31 //   I make a fit through the parameterization of the  FitRiemanTilt -> in order to get a straight line (representing the track)
  32 //   the distance of each cluster with respect to the straight line is calculated
  33 //   If there is no cluster at a distance of 0.6  to 2 with respect to the track, the track is defined as a good track, i.e. clean single track
  34 //                                                                        //
  35 //  Authors: Yvonne C. Pachmayer <pachmay@physi.uni-heidelberg.de>        //
  36 //                                                                        //
  37 ////////////////////////////////////////////////////////////////////////////
  38
  39 #include <TClonesArray.h>
  40 #include <TObjArray.h>
  41 #include <TProfile.h>
  42 #include <TH1.h>
  43 #include <TH1F.h>
  44 #include <TFile.h>
  45 #include "AliTRDtracker.h"
  46 #include "TTreeStream.h"
  47
  48 #include "AliPID.h"
  49 #include "AliESDtrack.h"
  50 #include "AliTrackReference.h"
  51 #include "AliExternalTrackParam.h"
  52 #include "AliTracker.h"
  53 #include "AliTRDseedV1.h"
  54 #include "AliTRDtrackV1.h"
  55 #include "AliTRDtrackerV1.h"
  56 #include "AliTrackPointArray.h"
  57 //#include "AliMagFMaps.h"
  58 #include "AliTRDcluster.h"
  59 #include "AliAnalysisManager.h"
  60
  61 #include "Cal/AliTRDCalPID.h"
  62 #include "AliTRDgeometry.h"
  63 #include "AliTRDmultiplicity.h"
  64 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  65
  66 ClassImp(AliTRDmultiplicity)
  67
  68 //____________________________________________________________________
  69 AliTRDmultiplicity::AliTRDmultiplicity()
  70     :AliTRDrecoTask("Multiplicity", "Barrel Tracking Multiplicity")
  71   ,fEventCounter(0)
  72 {
  73   //
  74   // Default constructor
  75     //
  76 }
  77
  78 //____________________________________________________________________
  79 AliTRDmultiplicity::~AliTRDmultiplicity()
  80 {
  81 }
  82
  83 //____________________________________________________________________
  84 void  AliTRDmultiplicity::CreateOutputObjects()
  85 {
  86   //
  87   // Create output objects
  88   //
  89
  90   OpenFile(0, "RECREATE");
  91
  92   TH1 *h = 0x0;
  93   fContainer = new TObjArray();
  94   for(Int_t is=0; is<AliTRDgeometry::kNsector; is++){
  95   fContainer->Add(h = new TH1F(Form("h_sector_%i", is), " ", 100,-10,10));
  96   }
  97 }
  98
  99 //____________________________________________________________________
 100 void AliTRDmultiplicity::Exec(Option_t *)
 101 {
 102   //
 103   // Do it
 104   //
 105
 106
 107
 108   ULong_t status;
 109   AliTRDtrackInfo     *track = 0x0;
 110   AliExternalTrackParam *esd = 0x0;
 111   AliTRDtrackV1 *TRDtrack = 0x0;
 112   Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
 113   Int_t stand_alone=1;
 114   fEventCounter++;
 115
 116
 117   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
 118     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 119     if(!track->HasESDtrack()) continue;
 120     status = track->GetStatus();
 121
 122     AliTrackPoint points[6];
 123     Float_t xyz[3];
 124     memset(xyz, 0, sizeof(Float_t) * 3);
 125     Float_t tracklet_x[6];
 126     Float_t tracklet_y[6];
 127     Float_t tracklet_z[6];
 128     for(Int_t a=0;a<6;a++)
 129     {
 130         xyz[0] = x_anode[a];
 131         points[a].SetXYZ(xyz);
 132         tracklet_x[a]=-1000;
 133         tracklet_y[a]=-1000;
 134         tracklet_z[a]=-1000;
 135     }
 136     Int_t det_tracklet=600;
 137
 138
 139
 140     // TRD standalone
 141     if(((status&AliESDtrack::kTRDout)>0) && !((status&AliESDtrack::kTRDin)>0))
 142     {
 143         TRDtrack = track->GetTrack();
 144
 145         if(TRDtrack)
 146         {
 147             AliTRDseedV1 *tracklet = 0x0;
 148             Int_t counter_tracklet=0;
 149             for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++)
 150             {
 151                 tracklet = TRDtrack->GetTracklet(itl);
 152                 if(tracklet)
 153                 {
 154                     if(tracklet->IsOK())
 155                     {
 156                         counter_tracklet++;
 157                         det_tracklet=tracklet->GetDetector();
 158                         //printf("%i %f %f \n",itl,tracklet->GetYfit(0),tracklet->GetZfit(0));
 159                         tracklet_x[itl]=tracklet->GetX0();
 160                         tracklet_y[itl]=tracklet->GetYfit(0);
 161                         tracklet_z[itl]=tracklet->GetZfit(0);
 162                     }
 163                 }
 164             }
 165             // this cut is needed otherwise the order of tracklets in the fit function is not correct
 166             if(counter_tracklet==AliTRDgeometry::kNlayer) AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(TRDtrack), 0x0, kTRUE, counter_tracklet, points);
 167             else continue;
 168
 169             if(fDebugLevel>=1){
 170               for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 171                 //printf("---------------------------------------- %i %i %f %f %f \n",counter_tracklet,il,points[il].GetX(),points[il].GetY(),points[il].GetZ());
 172                 Double_t point_x=points[il].GetX();
 173                 Double_t point_y=points[il].GetY();
 174                 Double_t point_z=points[il].GetZ();
 175
 176
 177                 (*fDebugStream)   << "TrackletsinTRD"
 178                     << "standalone=" << stand_alone
 179                     << "eventcounter=" << fEventCounter
 180                     << "layer="     << il
 181                     << "dettracklet=" << det_tracklet
 182                     << "xtracklet=" << tracklet_x[il]
 183                     << "xtrack="    << point_x
 184                     << "ytracklet=" << tracklet_y[il]
 185                     << "ytrack="    << point_y
 186                     << "ztracklet=" << tracklet_z[il]
 187                     << "ztrack="    << point_z
 188                     << "num_tracklets=" << counter_tracklet
 189                     << "\n";
 190
 191               }
 192             }
 193         }
 194     } // end TRD standalone
 195
 196
 197
 198     // TPC cluster selection for cosmic data
 199     if((track->GetTPCncls())<40) continue;
 200     // missing TPC propagation
 201     if(!(status&AliESDtrack::kTPCout)) continue;
 202     stand_alone=0;
 203
 204     esd = track->GetESDinfo()->GetOuterParam();
 205     // calculate sector - does not matter if track ends in TPC or TRD - sector number independent
 206     Double_t alpha;
 207     alpha=-100;
 208     Int_t fSector=100;
 209     if(esd){
 210         alpha=esd->GetAlpha();
 211         fSector = (Int_t)((alpha+ TMath::Pi())/(20*TMath::Pi()/180));
 212         if((fSector>-1) && (fSector<18))
 213         {
 214         }
 215         else
 216         {
 217             continue;
 218         }
 219     }
 220
 221
 222
 223     // only these sectors have a TRD detector at the moment
 224    if(fSector==0||fSector==8||fSector==9||fSector==17)
 225     {
 226         TRDtrack = track->GetTrack();
 227         if(!TRDtrack) continue;
 228         AliTRDseedV1 *tracklet = 0x0;
 229         Int_t counter_tracklet=0;
 230
 231
 232         for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 233             tracklet = TRDtrack->GetTracklet(itl);
 234             if(!tracklet || !tracklet->IsOK()) continue;
 235             counter_tracklet++;
 236             det_tracklet=tracklet->GetDetector();
 237             tracklet_x[itl]=tracklet->GetX0();
 238             tracklet_y[itl]=tracklet->GetYfit(0);
 239             tracklet_z[itl]=tracklet->GetZfit(0);
 240
 241             /*
 242             AliTRDcluster *c = 0x0;
 243             for(Int_t itime = 0; itime < AliTRDseedV1::kNtb; itime++){
 244                 c = tracklet->GetClusters(itime);
 245                 if(!c) continue;
 246 //                Float_t cluster_x = c->GetX();
 247 //                Bool_t isprop;
 248 //                isprop=kFALSE;
 249 //                isprop=AliTracker::PropagateTrackTo(esd,(Double_t)cluster_x,0.105,5,kFALSE);
 250 //                if(isprop)
 251 //                {
 252 //                    Int_t detector = c->GetDetector();
 253 //                    ((TH1F*)fContainer->At(fSector))->Fill((c->GetY())-(esd->GetY()));
 254 //                    if(c->GetY()-esd->GetY()>);
 255 //                    printf("diff: %f\n",c->GetY()-(esd->GetY()));
 256 //                }
 257             }
 258             */
 259
 260         } // loop over tracklets ends
 261
 262         if(counter_tracklet==AliTRDgeometry::kNlayer)
 263         {
 264             AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(TRDtrack), 0x0, kTRUE, counter_tracklet, points);
 265         }
 266         else continue;
 267
 268
 269
 270         if(fDebugLevel>=1){
 271           for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 272             //printf("---------------------------------------- %i %i %f %f %f \n",counter_tracklet,il,points[il].GetX(),points[il].GetY(),points[il].GetZ());
 273             Double_t point_x=points[il].GetX();
 274             Double_t point_y=points[il].GetY();
 275             Double_t point_z=points[il].GetZ();
 276
 277             (*fDebugStream)   << "TrackletsinTRD"
 278                 << "standalone=" << stand_alone
 279                 << "eventcounter=" << fEventCounter
 280                 << "layer="     << il
 281                 << "dettracklet=" << det_tracklet
 282                 << "xtracklet=" << tracklet_x[il]
 283                 << "xtrack="    << point_x
 284                 << "ytracklet=" << tracklet_y[il]
 285                 << "ytrack="    << point_y
 286                 << "ztracklet=" << tracklet_z[il]
 287                 << "ztrack="    << point_z
 288                 << "num_tracklets=" << counter_tracklet
 289                 << "\n";
 290           }
 291         }
 292
 293     }
 294   }
 295
 296
 297   PostData(0, fContainer);
 298 }
 299
 300 //____________________________________________________________________
 301 void AliTRDmultiplicity::Terminate(Option_t *)
 302 {
 303   //
 304   // Terminate
 305   //
 306
 307   if(fDebugStream){
 308     delete fDebugStream;
 309     fDebugStream = 0x0;
 310     fDebugLevel = 0;
 311   }
 312
 313   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 314   if (!fContainer) {
 315     Printf("ERROR: list not available");
 316     return;
 317   }
 318
 319 }
 320
 321
 322